/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Public Attributes | Static Public Attributes | Private Member Functions | Private Attributes
RollingGain Class Reference

#include <ClassRollingGain.h>

Collaboration diagram for RollingGain:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 RollingGain (string)
int AddResidual (string)
int AddNewHistogram (string)
int TestFit (double)
int Close ()
int LastFitProcess (string, string, string, bool)
void SetFitLast (bool)

Public Attributes

TimeStamp myBegin
TimeStamp myEnd

Static Public Attributes

static const string DetNames [14] = {"DayaBayAD1","DayaBayAD2","DayaBayIWS","DayaBayOWS","LingAoAD1","LingAoAD2","LingAoIWS","LingAoOWS","FarAD1","FarAD2","FarAD3","FarAD4","FarIWS","FarOWS"}

Private Member Functions

int Fit (DetProp &)
int Remove_hAdc (DetProp &)

Private Attributes

FILE * MasterFile
map< Detector, DetPropDetPropMap
TFile * ResidualDataFile
int low
int high
int nbins
const double PedestalWidth
TF1 * f1
TF1 * f2
string run1
string runf
bool FitLast

Detailed Description

Definition at line 28 of file ClassRollingGain.h.


Constructor & Destructor Documentation

RollingGain::RollingGain ( string  FileListName)

Definition at line 94 of file ClassRollingGain.h.

                                           :
        FitLast(false),
        PedestalWidth(3.8)
{

        f1 = new TF1("f1","gaus",20,250);
        f2 = new TF1("f2",NIMmodel,20,300,8);

        //set master file
        ifstream iMasterFile("Master.pmtCalibMap.txt");
        if(!iMasterFile){
                iMasterFile.close();
        MasterFile = fopen("Master.pmtCalibMap.txt", "w+");
        fprintf(MasterFile, "# [StartRun] [EndRun] [Sites] [SimFlag] [StartTimeSec] [StartTimeNanoSec] [EndTimeSec] [EndTimeNanoSec] [map file]\n");
        fclose(MasterFile);
        }
        else iMasterFile.close();

        //---------------Get Low and High --------------------------------
        low = 0;
        high = 0;
        int width = 900;
        ifstream FileList(FileListName.c_str());
        //First file, the beinning of histogram
        string NewHistogramFileName;
        getline(FileList,NewHistogramFileName);
        TFile* NewHistogramFile = new TFile( NewHistogramFileName.c_str(), "READ" );
        if( NewHistogramFile->IsOpen() ==0){
                cout << "ERROR: " << NewHistogramFileName << " doesn't exist!" << endl;
                exit(EXIT_FAILURE);     
        }
        run1 = NewHistogramFileName.substr(NewHistogramFileName.find("_run")+4,7);
        for( unsigned int idx=0; idx<14; idx++ )  {
      string BriefDetName = DetNames[idx];
      if( BriefDetName.find("_2In") != string::npos ) {
         BriefDetName = BriefDetName.substr( 0, BriefDetName.size()-4 );
      }
      string DetPathName = "detector_"+BriefDetName;
      TDirectory* pDetPath = (TDirectory*) NewHistogramFile->FindObjectAny( DetPathName.c_str() );
                if( !pDetPath )  continue;
                TObjLink *lnk = pDetPath->GetListOfKeys()->FirstLink();
        while( lnk ) {
                        string ChannelName = lnk->GetObject()->GetTitle();
        if( ChannelName.find("channel_") == string::npos ) {
                lnk = lnk->Next();
                continue;
        }
                string ChlPath = pDetPath->GetPath();
        ChlPath += "/"+ChannelName+"/dndadcVsTimeFine";
                TH2F* dndadcVsTimeFine = 0;
                dndadcVsTimeFine = (TH2F*) NewHistogramFile->Get( ChlPath.c_str() );
                           if( dndadcVsTimeFine ) {
                                int begin = (int) (dndadcVsTimeFine->GetXaxis()->GetXmin());
                                if( low ==0 ) low = begin;
                                else if( begin < low ) low = begin;
        }//end of if(dndadcVsTimeFine)
         lnk = lnk->Next();
                }//end of while(lnk)
        }
        NewHistogramFile->Close();
        //get the last
        while(!FileList.eof()){
                string temp;
                getline(FileList,temp);
                if(temp!="") NewHistogramFileName = temp;
   }
        //last file, the end of time
   NewHistogramFile = new TFile( NewHistogramFileName.c_str(), "READ" );
        runf = NewHistogramFileName.substr(NewHistogramFileName.find("_run")+4,7);
   for( unsigned int idx=0; idx<14; idx++ )  {
      string BriefDetName = DetNames[idx];
      if( BriefDetName.find("_2In") != string::npos ) {
         BriefDetName = BriefDetName.substr( 0, BriefDetName.size()-4 );
      }
      string DetPathName = "detector_"+BriefDetName;
      TDirectory* pDetPath = (TDirectory*) NewHistogramFile->FindObjectAny( DetPathName.c_str() );
                if( !pDetPath )  continue;
      TObjLink *lnk = pDetPath->GetListOfKeys()->FirstLink();
      while( lnk ) {
         string ChannelName = lnk->GetObject()->GetTitle();
         if( ChannelName.find("channel_") == string::npos ) {
            lnk = lnk->Next();
            continue;
         }
         string ChlPath = pDetPath->GetPath();
         ChlPath += "/"+ChannelName+"/dndadcVsTimeFine";
         TH2F* dndadcVsTimeFine = 0;
         dndadcVsTimeFine = (TH2F*) NewHistogramFile->Get( ChlPath.c_str() );
            if( dndadcVsTimeFine ) {
            int end = (int) (dndadcVsTimeFine->GetXaxis()->GetXmax());
            if( high == 0 ) high = end;
            high =  end > high ? end : high;
         }//end of if(dndadcVsTimeFine)
         lnk = lnk->Next();
      }//end of while(lnk)
   }
        low = low - 900;        //make sure all data are included.
   high = low + (int((high-0.1-low)/width)+1)*width;
   nbins = (high-low)/width; // 10 minutes;
        if(nbins>672) cout << "Warning: nbins indicates the data abnormally spread over 1 week!" << endl;
}

Member Function Documentation

int RollingGain::Fit ( DetProp currDet) [private]

Check which pmtDataTable to start by examing whether the file exists.

Pick out the detector need to fit

skip empty channels

Default values are ok to store.

skip low statisticd tubes

start a real fit

Step2: The gain fitting ^^^^^^^^^^^^^^^^^^^^^^^ Do a smarter fitting to automatically adjust some parameters until get a good fit So far it is found that lowerbound is the most sensitive parameter.

lowerbound trial is ordered by the chance a good fit can be reached.

uppperTrial is for the simple gauss fit only

Start right above pedestal then include the long tail.

Start right above pedestal then include the long tail.

Step3: Let's see if this channel can be bailed out by a simple gauss fit ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Start right above pedestal then include the long tail.

Last step: Fit result summary ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Definition at line 419 of file ClassRollingGain.h.

                                      {

        int NStop = 1;
        TFile * RootOutput = new TFile( "AdcFit.root", "UPDATE" );

        ostringstream *tmp = new ostringstream;
        *tmp <<setw(2) <<setfill('0') << NStop;
        string filename = "pmtDataTable_"+ tmp->str() +".txt";
        ifstream *ipmtTable = new ifstream(filename.c_str());
        while(ipmtTable->good()){
        ipmtTable = 0;
        NStop ++;
        tmp = new ostringstream; 
        *tmp <<setw(2) <<setfill('0') << NStop;
        filename = "pmtDataTable_"+tmp->str()+".txt";
        ipmtTable = new ifstream(filename.c_str());
        }       
        ipmtTable->close();  

    
        TimeStamp FirstTime = currDet.FitBegin;
        TimeStamp LastTime = currDet.FitEnd;
        
        FILE* m_outfile = fopen(filename.c_str(), "w+");
        cout << endl <<filename << " has been created." << endl;  

        fprintf(m_outfile, "# First Trigger time: (TrigSecond, TrigNanoSec) = (%d,%d).\n", (int)FirstTime.GetSec(), 
         (int)FirstTime.GetNanoSec());
        fprintf(m_outfile, "# Last Trigger time: (TrigSecond, TrigNanoSec) = (%d,%d).\n", (int)LastTime.GetSec(), 
          (int)LastTime.GetNanoSec()+1);
        fprintf(m_outfile, "# [pmtID] [channelID] [description] [status] [nhits] [fitStat] [chi2ndf] [speHigh] [gainErr] [sigmaSpe] [speLow] [timeOffset] [timeSpread] [efficiency] [prePulse] [afterPulse] [darkRate] [darkRateErr] [elecNoiseRate] [elecNoiseRateErr] [preAdc] [preAdcErr]\n");
  
        RootOutput->cd();
        string newdir = tmp->str();
        RootOutput->mkdir(newdir.c_str());
        RootOutput->cd(newdir.c_str());
  
        MasterFile = fopen("Master.pmtCalibMap.txt", "a");
        fprintf(MasterFile, "0000000 0000000 %4.0d %4.0d %10.0d 0 %10.0d 0 %s\n",
                         (currDet.Det.fullPackedData()>>16), 1, 
                         (int)FirstTime.GetSec(),
                         (int)LastTime.GetSec(),  filename.c_str());
        fclose(MasterFile);


        map<ElecChannelId,PmtProp>::iterator it,idend=currDet.PmtPropMap.end();
        for(it=currDet.PmtPropMap.begin(); it!=idend; it++) {
                PmtProp& curr = it->second;
        curr.Status=PmtCalibData::kUnknown;


                if(curr.h_Adc == 0){
                        curr.Status=PmtCalibData::kUnknown;

                }

        else if( curr.h_Adc->GetEntries() == 0 ) {          

        curr.Status=PmtCalibData::kUnknown;



        } 



                else if ( curr.h_Adc->GetEntries() > 0 &&    
                          curr.h_Adc->GetEntries() < 100 ) {

        curr.Status = PmtCalibData::kBadEfficiency;
        curr.NHits = (int)curr.h_Adc->GetEntries();
        //curr.FitStatus = 
        //curr.Chi2 = 
        curr.Gain = curr.h_Adc->GetMean();
        curr.Sigma = curr.h_Adc->GetRMS();
        curr.GainErr = curr.Sigma/TMath::Sqrt(curr.NHits);

        } 


                else {                                          
        ostringstream tmp3;
        tmp3<<setw(2)<<setfill('0')<< NStop;

        // Step1: fit pedestal
        // ^^^^^^^^^^^^^^^^^^^
        double q0 = 0;                  /* Pedestal subtracted. */
        double sigma0 = PedestalWidth;  /* I can't determine this here. */
        //----------------------------------------
        double AreaGuess, Q1MeanGuess, Q1SigmaGuess;
        AreaGuess = curr.h_Adc->GetEntries();
        Q1MeanGuess = curr.h_Adc->GetMean();
        Q1SigmaGuess = curr.h_Adc->GetRMS();
      
        double lowerbound, upperbound;
        double lowerTrial[9]={0.65*Q1MeanGuess,0.625*Q1MeanGuess,0.6*Q1MeanGuess,
                                                                         0.55*Q1MeanGuess,0.5*Q1MeanGuess,0.45*Q1MeanGuess,
                                                                         0.4*Q1MeanGuess,0.35*Q1MeanGuess,0.3*Q1MeanGuess};
      
        double upperTrial[9]={1.65*Q1MeanGuess,1.55*Q1MeanGuess,2.00*Q1MeanGuess,
                                 1.65*Q1MeanGuess,1.55*Q1MeanGuess,2.00*Q1MeanGuess,
                                 1.65*Q1MeanGuess,1.55*Q1MeanGuess,2.00*Q1MeanGuess};

        bool GotGoodFit = false;

        int status2;      // Fit status
        double q1;        // Gain
        double sigma1;    // Sigma
        double q1err;     // Gain Err
        double sigma1err; // Sigma Err
        double avechi2;   // Chi2/ndf

        for( unsigned int tidx=0; tidx<9; tidx++ ) {
                                lowerbound = q0+lowerTrial[tidx];    
                                upperbound = q0+Q1MeanGuess+5*Q1SigmaGuess;

                                // Initialize fitting parameters
                                f2->SetParameters(AreaGuess*10, 0, 0, Q1MeanGuess, Q1SigmaGuess, 0.01, 0.03, 0.1);
                                f2->SetParNames("N","Q_{0}","#sigma_{0}","Q_{1}","#sigma_{1}","w","#alpha","#mu");
      
                                f2->FixParameter(1,q0);
                                f2->FixParameter(2,sigma0);
                                f2->FixParameter(5,0.);
                                f2->FixParameter(6,0.);
                                //f2->FixParameter(5,w);
                                //f2->FixParameter(6,a);
      
                                f2->SetRange(lowerbound, upperbound); 
                                status2 = curr.h_Adc->Fit("f2","RQ");
        
                                q1 = f2->GetParameter(3);
                                sigma1 = f2->GetParameter(4);
                                q1err = f2->GetParError(3);
                                sigma1err = f2->GetParError(4);
                                avechi2 = f2->GetChisquare()/f2->GetNDF();
                                if(sigma1<0) sigma1 = -sigma1;//sometimes fit gives a negative sigma
        
                                // Keep trying until get a good fit.
                                // Define a good fit
                                // ^^^^^^^^^^^^^^^^^
                                if( status2==0 &&  // Root TH1 fit status 0 if ok (http://root.cern.ch/root/html/TH1.html#TH1:Fit)
                                avechi2<10 &&                    // Good Chi2/ndf
                                q1>6 && q1<80 &&                 // Too low or too high gain
                                q1>lowerbound && q1<upperbound &&// Within lower and upper bounds otherwise it always has big error
                                q1err<10 &&                     // Reasonable err of Gain
                                sigma1<q1 &&                     // Gain sigma is usually 1/3 of Gain
                                sigma1>1 && sigma1<30){          // Too low or too high sigma

                                        // Escape from the loop
                                        GotGoodFit = true;
                                        break;
                                }
                }




        if( !GotGoodFit ) {
                                cout<<"Batch "<<tmp3.str()<<" bad gain fit for PMT: "<<
                                curr.PmtLabel<<" try a simple Gaussian fit." <<endl;
        
                                for( unsigned int tidx=0; tidx<9; tidx++ ) {
                                        lowerbound = q0+lowerTrial[tidx];    
                                        upperbound = q0+upperTrial[tidx];

                                        f1->SetParameters(AreaGuess,Q1MeanGuess,7.);
                                        f1->SetRange(lowerbound, upperbound);
                                        status2 = curr.h_Adc->Fit("f1","RQ");
          
                                        q1 = f1->GetParameter(1)-q0;
                                        sigma1 = f1->GetParameter(2);
                                        q1err = f1->GetParError(1);
                                        sigma1err = f1->GetParError(2);
                                        avechi2 = f1->GetChisquare()/f1->GetNDF();
          
                                        // Define a good fit
                                        // ^^^^^^^^^^^^^^^^^
                                        if( status2==0 &&  // Root TH1 fit status 0 if ok (http://root.cern.ch/root/html/TH1.html#TH1:Fit)
                                avechi2<10 &&                         // Good Chi2/ndf
                                q1>6 && q1<80 &&                      // Too low or too high gain
                                q1>lowerbound && q1<upperbound && // Within lower and upper bounds otherwise it always has big error
                                q1err<10 &&                           // Reasonable err of Gain
                                sigma1<q1 &&                          // Gain sigma is usually 1/3 of Gain
                                sigma1>1 && sigma1<40                 // Too low or too high sigma
                        ) {
                                        // Escape from the loop
                                        GotGoodFit = true;
                                        break;
                                        }
                                }
        }




        if( !GotGoodFit ) {
                                cout << "Batch "<<tmp3.str()<<" gaussian fit for PMT: "
                                  << curr.PmtLabel
                                  <<" is not good either. Set the gain with mean value of h_Adc."<<endl;
                                q1 = curr.h_Adc->GetMean();
                                q1err = curr.h_Adc->GetMean() / TMath::Sqrt( curr.h_Adc->GetEntries() );
                                sigma1 = 0;
                                sigma1err = 0;
        }



        //q0 is  always 0 here
        //curr.Pedestal = q0; 
        curr.Gain = q1;
        curr.Sigma = sigma1;
        curr.GainErr = q1err;
        curr.SigmaErr = sigma1err;
        curr.NHits = (int)AreaGuess;
        curr.FitStatus = status2;
        curr.Chi2 = avechi2;
        if( !GotGoodFit ) {
                                curr.Status |= PmtCalibData::kBadGain;
        }
        if( curr.Sigma/curr.Gain > 1 ) {
                                curr.Status |= PmtCalibData::kBadGainWidth;
        }
        if( curr.Status == PmtCalibData::kUnknown) {
                                curr.Status = PmtCalibData::kGood;
        }
        }//else(non-empty, enough counts)
    
        curr.PmtId = 0;
                if(curr.h_Adc == 0) {
                        continue; //skip empty ones 
                }
        if(curr.h_Adc->GetEntries()>0){//skip 0s    


        fprintf(m_outfile, "%d %d %s %2d %6d %2d %6.2f %8.3f %8.3f  %7.3f %.1f %.1f %.1f %.1f %.1f %.1f %9.1f %9.1f %8.3f %8.3f %8.3f %8.3f\n", 
              curr.PmtId, curr.ChannelId, "none", curr.Status, curr.NHits, curr.FitStatus,
              curr.Chi2, curr.Gain, curr.GainErr, curr.Sigma, curr.Gain/19 ,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.);
      
        curr.h_Adc->Write();
        }//if(curr.h_Adc->GetEntries()>0)



        }//for(pmt)
   


  RootOutput->Close();
  fclose(m_outfile);

  return NStop;
}
int RollingGain::Remove_hAdc ( DetProp ADetector) [private]

Definition at line 384 of file ClassRollingGain.h.

                                              {
        for( map<ElecChannelId ,PmtProp>::iterator it = ADetector.PmtPropMap.begin();
                  it != ADetector.PmtPropMap.end(); it++ )  
        {
        PmtProp& APmt = it->second;
      delete APmt.h_Adc;
      APmt.h_Adc = 0 ;
   }
        return 0;
}
int RollingGain::AddResidual ( string  ResidualDataFileName)

Definition at line 197 of file ClassRollingGain.h.

                                                       {
        // Create a residual file. All dndadc will be saved here for backup.
        ResidualDataFile = new TFile( ResidualDataFileName.c_str(), "RECREATE" );
        cout << "Create residual data file: " << ResidualDataFileName << endl;
        for( unsigned int idx=0; idx<14; idx++ ) {
        DetProp& ADetector = DetPropMap[ DetNames[idx] ];
        ADetector.Name     = DetNames[idx];
        ADetector.Begin    = TimeStamp(low,0);
        ADetector.End      = TimeStamp(high,0);
        ADetector.Statistics = 0;
        ADetector.Det      = Detector::from_short( Detector::siteDetPackedFromString( ADetector.Name ) );
                ResidualDataFile->mkdir( DetNames[idx].c_str() );
        ResidualDataFile->cd();
        }//End of detector

        return 1;
}
int RollingGain::AddNewHistogram ( string  NewHistogramFileName)

Definition at line 216 of file ClassRollingGain.h.

                                                           {
        /* Open the new histogram file */
        /* 222222222222222222222222222 */
        // Channel Path looks like
        //     kup_v45_run0034263_seq0113.root:/stats/diagnostics/run_0034263/detector_DayaBayOWS/
        //               channel_board07_connector16/dndadcVsTimeFine

        cout <<"Loading "<<NewHistogramFileName<<endl;
        if( NewHistogramFileName.size() == 0 )  {
                cout << " New Histogram File is empty " << endl;
                return 0;
        }
   TFile* NewHistogramFile = new TFile( NewHistogramFileName.c_str(), "READ" );
        if( NewHistogramFile->IsOpen() ==0){
      cout << "ERROR: " << NewHistogramFileName << " doesn't exist!" << endl;
      exit(EXIT_FAILURE);  
   }
        myBegin = TimeStamp(0,0);       // To see each file's begin and end.
        myEnd = TimeStamp(0,0); 
   /* Load new histograms */
        for( unsigned int idx=0; idx<14; idx++ )  {
        /* Get a Detector Property */
        DetProp& ADetector = DetPropMap[ DetNames[idx] ];
        /* Brief detector name with 2In removed */
        string BriefDetName = DetNames[idx];
        if( BriefDetName.find("_2In") != string::npos ) {
                BriefDetName = BriefDetName.substr( 0, BriefDetName.size()-4 );
                }
        /* The dir for this detector exists? */
        string DetPathName = "detector_"+BriefDetName;
        TDirectory* pDetPath = (TDirectory*) NewHistogramFile->FindObjectAny( DetPathName.c_str() );
        if( !pDetPath )  continue;
        // Loop over each channel
        int begin = 0;  //begin and end time for this detector
                int end = 0;
        TObjLink *lnk = pDetPath->GetListOfKeys()->FirstLink();
                while( lnk ) {
                string ChannelName = lnk->GetObject()->GetTitle();
                if( ChannelName.find("channel_") == string::npos ) {
                                lnk = lnk->Next();
                                continue;
                        }
                        Detector siteDet;
         short int board,channel;
         siteDet = Detector::from_short( Detector::siteDetPackedFromString( BriefDetName ) );
         board   = atoi( ChannelName.substr( ChannelName.find("board")+5, 2).c_str() );
         channel = atoi( ChannelName.substr( ChannelName.find("connector")+9, 2).c_str() );
         // This channel's ElecChannelId
         ElecChannelId ChannelId( board, channel, siteDet.site(), siteDet.detectorId() );
                // Its histogram
                string ChlPath = pDetPath->GetPath();
                ChlPath += "/"+ChannelName+"/dndadcVsTimeFine";
                // Get the right histogram
                TH2F* dndadcVsTimeFine = 0;
                dndadcVsTimeFine = (TH2F*) NewHistogramFile->Get( ChlPath.c_str() );
                if( dndadcVsTimeFine ) {
            string HistoName = DetNames[idx]+"_"+ChannelName;
                                //find begin and end time of this detector
                int tempBegin = (int) (dndadcVsTimeFine->GetXaxis()->GetXmin());
                int tempEnd   = (int) (dndadcVsTimeFine->GetXaxis()->GetXmax());
                                if(begin==0) begin = tempBegin; end = tempEnd;
                                begin = tempBegin < begin? tempBegin : begin;
            end = tempEnd > end? tempEnd : end;
                                //boundary check
                                if(tempBegin<low || tempEnd>high) cout<< "Warning: "<<ChlPath<<" extends out of boundary! "<<endl;
                                //Prepare Residual File
                                PmtProp &APmt = ADetector.PmtPropMap[ ChannelId ];
                                if( !APmt.h_AdcVsTime ){
                                        //Residual file doesn't exist. Create one.
                                        APmt.h_AdcVsTime = new TH2F(HistoName.c_str(),HistoName.c_str(),nbins,low,high,480,-5,115);
                                        APmt.ChannelId = ChannelId.fullPackedData();
                                }
                                //Add new histogram to residual data
                                for(int binx = 1;binx <= dndadcVsTimeFine->GetXaxis()->GetNbins();binx++ )  {
                                for(int biny = 1;biny <= dndadcVsTimeFine->GetYaxis()->GetNbins();biny++ )  {
                                double x = dndadcVsTimeFine->GetXaxis()->GetBinCenter( binx );
                                double y = dndadcVsTimeFine->GetYaxis()->GetBinCenter( biny );
                                      int inc = (int) dndadcVsTimeFine->GetBinContent( binx,biny );
                                for(int i = 0 ; i < inc ; i ++){
                                        APmt.h_AdcVsTime->Fill(x,y);
                                }
                                }
                                }
            ADetector.Statistics += (int) dndadcVsTimeFine->Integral();
                                APmt.h_AdcVsTime->SetDirectory(0);
         }//end of if(dndadcVsTimeFine)
         lnk = lnk->Next();

      }//end of while(lnk)
                
                //Get the begin and end time for this file.
      //ADetector.Begin = TimeStamp(begin,0);
      ADetector.End   = TimeStamp(end,0);
                if(myBegin.GetSec()==0) myBegin = TimeStamp(begin,0) ; myEnd = TimeStamp(end,0);
                myBegin = begin < myBegin.GetSec() ? TimeStamp(begin,0) : myBegin;
                myEnd = end > myEnd.GetSec() ? TimeStamp(end,0) : myEnd;                
        }//end of for(dx)
        NewHistogramFile->Close();
        return 1;
}
int RollingGain::TestFit ( double  FitStatistic)

Definition at line 318 of file ClassRollingGain.h.

                                           {
   for( unsigned int idx=0; idx<14; idx++ )  {
      /* Get a Detector Property */
      DetProp& ADetector = DetPropMap[ DetNames[idx] ];
      int nChannels = ADetector.PmtPropMap.size();
      if( nChannels == 0) continue;
      double aveStat = (nChannels == 0) ? 0 : ADetector.Statistics/nChannels;
      int Start = 1;
      bool lastfit = FitLast;
                //output aveStat situation
      if(nChannels!= 0) {
         cout << DetNames[idx] << setw(15) <<  aveStat << setw(15) << "From " << setw(15) << ADetector.Begin 
                                  << setw(15) << " To " << setw(15) << ADetector.End  << endl;
      }
      // This det has more events than expected.
      while( aveStat > FitStatistic || lastfit ) {
         // Start to make projections since the beginning.
         for( int BinIdx = Start; BinIdx <= nbins; BinIdx++ )  {
            // Make Projections for all channels for this detector
            int CurrentStat = 0;
            for( map<ElecChannelId ,PmtProp>::iterator it = ADetector.PmtPropMap.begin();
                 it != ADetector.PmtPropMap.end(); it++ )  {
               PmtProp& APmt = it->second;
               APmt.h_Adc = APmt.h_AdcVsTime->ProjectionY( "", Start, BinIdx );
               string ChlName = APmt.h_AdcVsTime->GetName();
               APmt.h_Adc->SetName( (ChlName+"_Adc").c_str() );
               APmt.h_Adc->SetTitle( (ChlName+"_Adc").c_str() );
               CurrentStat += (int) APmt.h_Adc->Integral();
               ADetector.FitBegin = TimeStamp( time_t( APmt.h_AdcVsTime->GetXaxis()->GetBinLowEdge(Start) ), 0);
               ADetector.FitEnd   = TimeStamp( time_t( APmt.h_AdcVsTime->GetXaxis()->GetBinUpEdge(BinIdx) ), 0);
            }
            double thisstat = (nChannels == 0) ? 0 : CurrentStat/nChannels;
                                // End of histogram
                                if( BinIdx == nbins )  {
               lastfit = false;
               cout<<"Last fit"<<endl;
               if(CurrentStat!=0) {
                                                Fit( ADetector );
                                                ADetector.Statistics = 0;
                                        }
               else {
                                                cout << "Last fit of " << DetNames[idx] << "has 0 counts. Fitting will not proceed." << endl;
                                                ADetector.Statistics = 0;
                                        }
               Remove_hAdc( ADetector );
            }
                                //Or enough fit
                                else if( thisstat > FitStatistic ){     
               Fit( ADetector );
               ADetector.Statistics -= CurrentStat;
               Start = BinIdx+1;
               Remove_hAdc( ADetector );
               break;
                                }
                                //Nothing
                                else Remove_hAdc( ADetector );
         }
         aveStat = (nChannels == 0) ? 0 : ADetector.Statistics/nChannels;
      }//while(aveStat > FitStatistic || lastfit)
      lastfit = FitLast;
   }//end of for(idx)

   return 1;
}
int RollingGain::Close ( )

Definition at line 395 of file ClassRollingGain.h.

                      {
//Close will save h_AdcVsTime to Residual file.
        for( unsigned int idx=0; idx<14; idx++ )  {
                TList *SaveList = new TList();
           /* Get a Detector Property */
        DetProp& ADetector = DetPropMap[ DetNames[idx] ];
        ResidualDataFile->cd();
        ResidualDataFile->cd( DetNames[idx].c_str() );
        for( map<ElecChannelId ,PmtProp>::iterator it = ADetector.PmtPropMap.begin();
                it != ADetector.PmtPropMap.end(); it++ )  {
        PmtProp& APmt = it->second;
        // TObject::kOverwrite will force a replacement of the previous object with the same name.
                        SaveList->Add(APmt.h_AdcVsTime);
                }
                SaveList->Write("",TObject::kOverwrite);        
        }
        ResidualDataFile->Close();
        return 1;
}
int RollingGain::LastFitProcess ( string  site,
string  runMode,
string  db,
bool  isqsubFile 
)

Definition at line 60 of file ClassRollingGain.h.

                                                                                      {
        /*
        1. Make a folder name run1-run2-time1-time2.
        2. mv pmtDataTable, MasterTable, AdcFit.root, Timer.txt
        */
        TimeStamp ti(low,0);
        TimeStamp tf(high,0);
        stringstream ss_ti;
        stringstream ss_tf;
        ss_ti << setw(4) << setfill('0') << ti.GetDate()%10000;
        ss_tf << setw(4) << setfill('0') << tf.GetDate()%10000;

        string dir = run1 + "-" + runf + "-" + ss_ti.str() + "-" + ss_tf.str();
   string mkdir = "mkdir " + dir;
   string mv = "mv *.root *.txt qsub.* ";
   mv = mv +  " ./" + dir;
   system(mkdir.c_str());
   system(mv.c_str());
        string nuwaCommand = "nuwa.py --dbconf=" + db + " -m'RollingGainAuto.ScanFrames " + runMode + " False'";
        string readout;
        readout = exec( nuwaCommand.c_str() );
        cout << endl << endl << readout << endl;
        system("rm frames.*");
        if(!isqsubFile){
                mv = "mv " + dir + " ../pass-" + site + "/";
                system( mv.c_str() );
        }
        return 1;
}
void RollingGain::SetFitLast ( bool  fitstatus)

Definition at line 92 of file ClassRollingGain.h.

{FitLast = fitstatus;}

Member Data Documentation

FILE* RollingGain::MasterFile [private]

Definition at line 29 of file ClassRollingGain.h.

Definition at line 30 of file ClassRollingGain.h.

Definition at line 31 of file ClassRollingGain.h.

int RollingGain::low [private]

Definition at line 34 of file ClassRollingGain.h.

int RollingGain::high [private]

Definition at line 35 of file ClassRollingGain.h.

int RollingGain::nbins [private]

Definition at line 36 of file ClassRollingGain.h.

const double RollingGain::PedestalWidth [private]

Definition at line 37 of file ClassRollingGain.h.

TF1* RollingGain::f1 [private]

Definition at line 38 of file ClassRollingGain.h.

TF1* RollingGain::f2 [private]

Definition at line 39 of file ClassRollingGain.h.

string RollingGain::run1 [private]

Definition at line 40 of file ClassRollingGain.h.

string RollingGain::runf [private]

Definition at line 41 of file ClassRollingGain.h.

bool RollingGain::FitLast [private]

Definition at line 42 of file ClassRollingGain.h.

Definition at line 50 of file ClassRollingGain.h.

Definition at line 51 of file ClassRollingGain.h.

const string RollingGain::DetNames = {"DayaBayAD1","DayaBayAD2","DayaBayIWS","DayaBayOWS","LingAoAD1","LingAoAD2","LingAoIWS","LingAoOWS","FarAD1","FarAD2","FarAD3","FarAD4","FarIWS","FarOWS"} [static]

Definition at line 52 of file ClassRollingGain.h.


The documentation for this class was generated from the following file:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 09:59:50 for RollingGainAuto by doxygen 1.7.4