loon: Trimmer Class Reference

00040 {
00041   
00042   NueRecord *nr = new NueRecord();
00043   NuePOT *np = new NuePOT();
00044 
00045   TChain *chain = new TChain("ana_nue");
00046   chain->SetBranchAddress("NueRecord",&nr);
00047                                                                         
00048   TChain *pchain = new TChain("pottree");
00049   pchain->SetBranchAddress("NuePOT", &np);
00050  
00051   for(unsigned int i = 0; i < files.size(); i++)
00052   {
00053       chain->Add(files[i].c_str());
00054       pchain->Add(files[i].c_str());
00055   }
00056   
00057   if(!outSet){
00058      string file;
00059                                                                         
00060      if(pchain->GetEntries() > 0){
00061         pchain->GetEntry(0);
00062         file = pchain->GetFile()->GetName();
00063      }
00064      string minifile = file.substr(file.find_last_of("/")+1, file.find_last_of(".root")-file.find_last_of("/") - 5);
00065      minifile += "-Trim.root";
00066   
00067      outFile = minifile;
00068   }
00069   if(outFile == "-Trim.root"){
00070      cout<<"No input file found"<<endl; 
00071      return;
00072   }
00073   
00074   string out1,out2,out3;
00075   if(separatebyRunPeriod)
00076   {
00077     string base = outFile.substr(0,outFile.size()-5);//strip off '.root' and add label for run period
00078     out1 = base + "_RunPeriod1.root";
00079     out2 = base + "_RunPeriod2.root";
00080     out3 = base + "_RunPeriod3.root";
00081     cout<<"Setting output to:"<<endl;
00082     cout<<out1<<endl;
00083     cout<<out2<<endl;
00084     cout<<out3<<endl;
00085   }
00086   else
00087   {
00088     cout<<"Setting output to "<<outFile<<endl;
00089     out1 = outFile;
00090   }
00091 
00092   // And now the actual looping
00093   
00094   TFile *f1 = new TFile(out1.c_str(),"RECREATE");
00095   TFile *f2=0,*f3=0;
00096   if(separatebyRunPeriod)
00097   {
00098     f2 = new TFile(out2.c_str(),"RECREATE");
00099     f3 = new TFile(out3.c_str(),"RECREATE");
00100   }
00101   
00102   f1->cd();
00103   TTree *tree1 = new TTree("ana_nue","ana_nue");
00104   TTree::SetBranchStyle(1);
00105   TBranch* br1 = tree1->Branch("NueRecord","NueRecord", &nr );
00106   br1->SetAutoDelete(kFALSE);
00107   
00108   TTree *tree2=0;
00109   TTree *tree3=0;
00110   TBranch* br2=0;
00111   TBranch* br3=0;
00112   if(separatebyRunPeriod)
00113   {
00114     f2->cd();
00115     tree2 = new TTree("ana_nue","ana_nue");
00116     br2 = tree2->Branch("NueRecord","NueRecord", &nr );
00117     br2->SetAutoDelete(kFALSE);
00118     
00119     f3->cd();
00120     tree3 = new TTree("ana_nue","ana_nue");
00121     br3 = tree3->Branch("NueRecord","NueRecord", &nr );
00122     br3->SetAutoDelete(kFALSE);
00123   }
00124   
00125   bool isRunPeriod1=false;
00126   bool isRunPeriod2=false;
00127   bool isRunPeriod3=false;
00128   
00129   Int_t n = chain->GetEntries();
00130   int count = 0;
00131 
00132   bool isMC = false;
00133   
00134   //for recalculating the pot tree for each run period
00135   double TotPOT1 = 0.0,TotPOT2 = 0.0,TotPOT3 = 0.0;
00136   double TotPOT1_nocut = 0.0,TotPOT2_nocut = 0.0,TotPOT3_nocut = 0.0;
00137   Int_t nruns=0;
00138   Int_t nsnarls1=0,nsnarls2=0,nsnarls3=0;
00139   
00140   int lastrun=-1;
00141   
00142   chain->GetEntry(0);
00143   ReleaseType::Release_t rel = nr->GetHeader().GetRelease();
00144   BeamType::BeamType_t beamtype = nr->GetHeader().GetBeamType();
00145   Detector::Detector_t det = nr->GetHeader().GetVldContext().GetDetector();
00146   
00147   for(int i=0;i<n;i++){
00148     if(i%10000==0) cout << 100*float(i)/float(n)
00149                        << "% done" << endl;
00150     chain->GetEvent(i);
00151     
00152     NueConvention::NueEnergyCorrection(nr);
00153     NueStandard::ModifyANNPID(nr);
00154     
00155     nr->eventq.passFarDetQuality = DataUtil::IsGoodData(nr->GetHeader().GetVldContext());
00156     nr->eventq.passNearDetQuality = nr->eventq.passFarDetQuality;
00157     
00158     if(i == 0){
00159       SimFlag::SimFlag_t sim = nr->GetHeader().GetVldContext().GetSimFlag();
00160       isMC = (sim == SimFlag::kMC);
00161     }
00162     
00163     if(separatebyRunPeriod && !isMC)
00164     {
00165       cout<<"'separatebyRunPeriod' option is not meant to be used for data - pottree would not be calculated correctly for data.  Qutting..."<<endl;
00166       return;
00167     }
00168     
00169     //sanity check
00170     if(NueStandard::PassesPOTStandards(nr))
00171     {
00172       if(nr->GetHeader().GetEventNo() == 0 || nr->GetHeader().GetEvents() == 0)
00173       {
00174         if(nr->bmon.bI < 0 && !isMC) cout<<"Unexpected behavior"<<endl;
00175       }
00176     }
00177    
00178     isRunPeriod1=false;
00179     isRunPeriod2=false;
00180     isRunPeriod3=false;
00181     if(separatebyRunPeriod)
00182     {
00183       if(nr->GetHeader().GetVldContext().GetTimeStamp().GetSec()<1145036420)
00184       {
00185         isRunPeriod1=true;
00186       }
00187       else if(nr->GetHeader().GetVldContext().GetTimeStamp().GetSec()<1190028111)
00188       {
00189         isRunPeriod2=true;
00190       }
00191       else
00192       {
00193         isRunPeriod3=true;
00194       }
00195     }
00196     
00197     if(separatebyRunPeriod)//recalculating pottree; this part works only for MC
00198     {
00199       if(nr->GetHeader().GetEventNo() == 0 || nr->GetHeader().GetEvents() == 0)//for each snarl
00200       {
00201         if(isRunPeriod1)
00202         {
00203           if(det==Detector::kNear)//this only works for ND MC
00204           {
00205             TotPOT1_nocut += MCInfo::GetMCPoT(det, beamtype, rel);
00206             TotPOT1 += MCInfo::GetMCPoT(det, beamtype, rel);
00207           }
00208           
00209           nsnarls1++;//count the snarls per run period
00210         }
00211         if(isRunPeriod2)
00212         {
00213           if(det==Detector::kNear)//this only works for ND MC
00214           {
00215             TotPOT2_nocut += MCInfo::GetMCPoT(det, beamtype, rel);
00216             TotPOT2 += MCInfo::GetMCPoT(det, beamtype, rel);
00217           }
00218           
00219           nsnarls2++;
00220         }
00221         if(isRunPeriod3)
00222         {
00223           if(det==Detector::kNear)//this only works for ND MC
00224           {
00225             TotPOT3_nocut += MCInfo::GetMCPoT(det, beamtype, rel);
00226             TotPOT3 += MCInfo::GetMCPoT(det, beamtype, rel);
00227           }
00228           
00229           nsnarls3++;
00230         }
00231       }
00232       
00233       if(nr->GetHeader().GetRun()!=lastrun)//counting runs - BUT not counting separately for each run period!!
00234       {
00235         lastrun = nr->GetHeader().GetRun();
00236         nruns++;
00237         
00238         if(det==Detector::kFar)//in FD MC, each run is a single file, and GetMCPoT returns pot/file
00239         {
00240           if(isRunPeriod1)
00241           {
00242             TotPOT1_nocut += MCInfo::GetMCPoT(det, beamtype, rel);
00243             TotPOT1 += MCInfo::GetMCPoT(det, beamtype, rel);
00244           }
00245           if(isRunPeriod2)
00246           {
00247             TotPOT2_nocut += MCInfo::GetMCPoT(det, beamtype, rel);
00248             TotPOT2 += MCInfo::GetMCPoT(det, beamtype, rel);
00249           }
00250           if(isRunPeriod3)
00251           {
00252             TotPOT3_nocut += MCInfo::GetMCPoT(det, beamtype, rel);
00253             TotPOT3 += MCInfo::GetMCPoT(det, beamtype, rel);
00254           }
00255         }
00256       }
00257     }
00258     
00259     if(fOverWritePOT && !isMC)//recalculating POT for data
00260     {
00261       if(NueStandard::PassesPOTStandards(nr))
00262       {
00263         if(nr->GetHeader().GetEventNo() == 0 || nr->GetHeader().GetEvents() == 0)
00264         {
00265           TotPOT1 += nr->bmon.bI;
00266         }
00267       }
00268     }
00269     
00270     if(EvaluateCuts(nr)){
00271        if(fReweight && nr->mctrue.nuFlavor > -10)
00272        {
00273           int nuFlavor = nr->mctrue.nuFlavor;
00274           int  nonOsc = nr->mctrue.nonOscNuFlavor;
00275           float energy = nr->mctrue.nuEnergy;
00276                                                                                 
00277           Float_t newWeight = NueConvention::Oscillate(nuFlavor, nonOsc, energy,                                735, fDeltaMSquare, fTheta23, fUe3Square);
00278                                                                                 
00279           nr->mctrue.Ue3Squared = fUe3Square;
00280           nr->mctrue.DeltamSquared23 = fDeltaMSquare;
00281           nr->mctrue.Theta23  = fTheta23;
00282           nr->mctrue.fOscProb = newWeight;
00283        }
00284        
00285        if(separatebyRunPeriod)
00286        {
00287          if(isRunPeriod1)
00288          {
00289            tree1->Fill();
00290          }
00291          if(isRunPeriod2)
00292          {
00293            tree2->Fill();
00294          }
00295          if(isRunPeriod3)
00296          {
00297            tree3->Fill();
00298          }
00299        }
00300        else
00301        {
00302          tree1->Fill();
00303        }
00304        count++;
00305     }
00306   }
00307   cout<<count<<" of "<<n<<" entries were passed"<<endl;
00308   
00309   f1->cd();
00310   NuePOT *total1 = new NuePOT();
00311   TTree *ptree1 = new TTree("pottree","pottree");
00312   TBranch* brp1 = ptree1->Branch("NuePOT","NuePOT", &total1 );
00313   brp1->SetAutoDelete(kFALSE);
00314   
00315   NuePOT *total2=0,*total3=0;
00316   TTree *ptree2=0,*ptree3=0;
00317   TBranch *brp2=0,*brp3=0;
00318   if(separatebyRunPeriod)
00319   {
00320     f2->cd();
00321     total2 = new NuePOT();
00322     ptree2 = new TTree("pottree","pottree");
00323     brp2 = ptree2->Branch("NuePOT","NuePOT", &total2 );
00324     brp2->SetAutoDelete(kFALSE);
00325     
00326     f3->cd();
00327     total3 = new NuePOT();
00328     ptree3 = new TTree("pottree","pottree");
00329     brp3 = ptree3->Branch("NuePOT","NuePOT", &total3 );
00330     brp3->SetAutoDelete(kFALSE);
00331   }
00332   
00333   if(pchain->GetEntries() > 0){
00334     pchain->GetEntry(0);
00335     total1->beamtype = np->beamtype;
00336     if(separatebyRunPeriod)
00337     {
00338       total2->beamtype = np->beamtype;
00339       total3->beamtype = np->beamtype;
00340     }
00341   }
00342   
00343   for(int i = 0; i < pchain->GetEntries(); i++)
00344   {
00345      pchain->GetEntry(i);
00346      
00347      if(total1->beamtype != np->beamtype)
00348      {
00349        cerr<<"You are merging files of different beamtype - BAD"<<endl;
00350      }
00351      
00352      total1->nruns +=  np->nruns;
00353      total1->nsnarls +=  np->nsnarls;
00354      total1->pot += np->pot;
00355      total1->pot_nocut += np->pot_nocut;
00356   }
00357 
00358   if(fOverWritePOT && !isMC) total1->pot = TotPOT1;
00359   
00360   if(separatebyRunPeriod)//recalculating pottree
00361   {
00362     total1->pot = TotPOT1;
00363     total2->pot = TotPOT2;
00364     total3->pot = TotPOT3;
00365     
00366     total1->pot_nocut = TotPOT1_nocut;
00367     total2->pot_nocut = TotPOT2_nocut;
00368     total3->pot_nocut = TotPOT3_nocut;
00369     
00370     total1->nruns = nruns;
00371     total2->nruns = nruns;
00372     total3->nruns = nruns;
00373     
00374     total1->nsnarls = nsnarls1;
00375     total2->nsnarls = nsnarls2;
00376     total3->nsnarls = nsnarls3;
00377   }
00378   
00379   //fill the pottree
00380   ptree1->Fill();
00381   if(separatebyRunPeriod)
00382   {
00383     ptree2->Fill();
00384     ptree3->Fill();
00385   }
00386   
00387   //save the file(s)
00388   f1->cd();
00389   tree1->Write("ana_nue",TObject::kWriteDelete);
00390   ptree1->Write("pottree",TObject::kWriteDelete);
00391   if(separatebyRunPeriod)
00392   {
00393     f2->cd();
00394     tree2->Write("ana_nue",TObject::kWriteDelete);
00395     ptree2->Write("pottree",TObject::kWriteDelete);
00396     
00397     f3->cd();
00398     tree3->Write("ana_nue",TObject::kWriteDelete);
00399     ptree3->Write("pottree",TObject::kWriteDelete);
00400   }
00401   
00402   f1->Close();
00403   if(separatebyRunPeriod)
00404   {
00405     f2->Close();
00406     f3->Close();
00407   }
00408   
00409   
00410   cout<<"Trimming completed with "<<np->pot<<"x10^12 POT"<<endl;
00411 }

Public Member Functions
	Trimmer ()
void	AddFiles (string files)
void	SetOutputFile (string name)
void	Config (Registry &r)
void	SetCuts (string type, int level=0)
bool	EvaluateCuts (NueRecord *nr)
void	RunTrimmer ()
void	SetDeltaMSquare (double dm2)
void	SetUe3Square (double dUe32)
void	SetTheta23 (double t23)
void	SetBaseline (double bl)
void	RecalculatePOT ()
void	SeparatebyRunPeriod ()
Private Attributes
string	outFile
bool	outSet
vector< string >	files
NueAnalysisCuts	fCuts
double	fDeltaMSquare
double	fUe3Square
double	fTheta23
double	fBaseline
bool	fReweight
string	cutSet
int	cutLevel
bool	fOverWritePOT
bool	separatebyRunPeriod
Trimmer Class Reference

Public Member Functions

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation
