Changeset 467 in svn for trunk/Resolutions.cpp

Timestamp:

Jul 10, 2009, 5:57:55 PM (15 years ago)

Author:

Xavier Rouby

Message:

cleaning, comments, progressbar, cout

File:

• trunk/Resolutions.cpp (modified) (13 diffs)

trunk/Resolutions.cpp

@@ -38 +38 @@
 #include "TApplication.h"
 #include "TFile.h"
+#include "TStopwatch.h"
 
 #include "ExRootTreeReader.h"
 #include "ExRootTreeWriter.h"
 #include "ExRootTreeBranch.h"
+#include "ExRootProgressBar.h"
 #include "TreeClasses.h"
 
@@ -189 +191 @@
 {
   int appargc = 2;
-  char *appName = "Resolution";
-  char *appargv[] = {appName, "-b"};
+  char *appName= new char[20];
+  char *appOpt= new char[20];
+  sprintf(appName,"Resolution_ATLAS");
+  sprintf(appOpt,"-b");
+  char *appargv[] = {appName,appOpt};
   TApplication app(appName, &appargc, appargv);
-  
+  delete [] appName;
+  delete [] appOpt;
+
   if(argc != 3) {
       cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
-      cout << " input_list - list of files in root format," << endl;
+      cout << " input_file - input file in Delphes-root format," << endl;
       cout << " output_file - output file." << endl;
       exit(1);
   }
 
-  srand (time (NULL));         /* Initialisation du générateur */
-  
-  //read the input TROOT file
+
+  // 1. ********** initialisation ***********
+
+  srand (time (NULL));
+  TStopwatch globalwatch;
+  globalwatch.Start();
+  
   string inputfilename(argv[1]), outputfilename(argv[2]);
 
+  // 1.1 checks the input file
+  if(inputfilename.find(".root") > inputfilename.length()   ) {
+    cout << "input_file should be a .root file from Delphes!\n";
+    return -1;
+  }
+  TFile f(inputfilename.c_str());
+  if (!f.FindKey("GEN") || !f.FindKey("Analysis")           ) {
+    cout << "input_file should be a .root file from Delphes!\n";
+    cout << "it should contain both \"GEN\" and \"Analysis\" trees at least.\n";
+    return -1;
+  }
+  f.Close();
+
+  // 1.2 checks the output file
   if(outputfilename.find(".root") > outputfilename.length() ) {
     cout << "output_file should be a .root file!\n";
     return -1;
   }
-
-
-
   TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
-  
+ 
+  // 1.3 Reads the trees in input file 
   TChain chainGEN("GEN");
   chainGEN.Add(inputfilename.c_str());
@@ -229 +252 @@
   TIter itGen((TCollection*)branchGen);
   
-  //write the output root file
+
+  // 1.4 Prepares the output root file
   ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
   ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
@@ -236 +260 @@
   ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
   ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
-  
+ 
+  // 1.5 other initialisations 
   TRootGenParticle *particle;
-  
   RESOLELEC * elementElec;
   RESOLMUON *elementMuon;
@@ -245 +269 @@
   ETMIS *elementEtmis;
   
-  int numTau=0;
-  int numTauRec=0;
+  int numTau=0, numTauRec=0;
   
   RESOLution *DET = new RESOLution();
+  /*
+  string detectorcard = "data/DetectorCard_CMS.dat";
+  const float dR_jetpairing = 0.25;
+  const float jet_pt_cut = 1;
+  */
   DET->ReadDataCard("data/DetectorCard_CMS.dat");
 
 
-  //Jet information
   JetsUtil *JETRUN = new JetsUtil("data/DetectorCard_CMS.dat");
 
-  TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
-  TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information
+  TLorentzVector genMomentum(0,0,0,0); //generator level information
+  TLorentzVector recoMomentum(0,0,0,0);//reconstructed level information
   LorentzVector jetMomentum;
   
   vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
   vector<fastjet::PseudoJet> sorted_jetsGEN;
- 
   vector<int> NTrackJet;
- 
   vector<TLorentzVector> towers;
  
-  // Loop over all events
+  // 2. Loop over all events
   Long64_t entry, allEntries = treeReader->GetEntries();
-  cout << "** Chain contains " << allEntries << " events" << endl;
+  cout << endl << endl;
+  cout <<"*********************************************************************"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**        ####### Start resolution processing    ########          **"<< endl;
+  cout << left  << setw(52) <<"**              Total number of events to run: "<<""
+       << left  << setw(15) << allEntries <<""
+       << right << setw(2) <<"**"<<endl;
+
+  ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
+
   for(entry = 0; entry < allEntries; ++entry)
     {
+      Progress->Update(entry);
       TLorentzVector PTmisReco(0,0,0,0);
       TLorentzVector PTmisGEN(0,0,0,0);
@@ -276 +313 @@
       treeReaderGEN->ReadEntry(entry);
       treeWriter->Clear();
-      if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
       
       TSimpleArray<TRootGenParticle> bGen;
@@ -301 +337 @@
           //Calculate ETMIS from generated particles
           if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
-         
+
+          //Electrons and muons  
           if( (particle->Status == 1)    &&
             ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
@@ -320 +357 @@
                  case pMU: // all muons with eta < DET->MAX_MU
                     PairingMuon(recoMomentum,genMomentum,branchMuon);
-                    if(recoMomentum.E() !=0){
+                    if(recoMomentum.E()!=0){
                      elementMuon = (RESOLMUON*) branchmuon->NewEntry();
-                     elementMuon->OverPT = (1/genMomentum.Pt());
-                     elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
+                     elementMuon->OverPT = 1./genMomentum.Pt();
+                     elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
                    break; // case pMU
                  default:
@@ -354 +391 @@
       elementEtmis->Ex = (-PTmisGEN).Px();
       elementEtmis->SEt = ScalarEt;
-
       elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
       elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
@@ -362 +398 @@
 
       TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
-
       TLorentzVector JETreco(0,0,0,0);
       for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
@@ -372 +407 @@
             elementJet= (RESOLJET*) branchjet->NewEntry();
             elementJet->PT = JETgen.Et();
-            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
+            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();  //// la difference pourrait être ici
           }
       }
-      numTau = numTau+TausHadr.GetEntries();
+      numTau += TausHadr.GetEntries();
 
       TIter itJet((TCollection*)branchJet);
@@ -391 +426 @@
                   {
                     elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
-                    elementTaujet->EnergieCen = (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E());
+                    elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
                     elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
                     if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95 
@@ -401 +436 @@
         
       } // for itJet : loop on all jets
-//cout<<"i"<<endl;
       
       treeWriter->Fill();
     } // Loop over all events
   treeWriter->Write();
-float frac = numTauRec/numTau;  
-cout<<numTauRec<<endl;
-cout<<numTau<<endl;
-
-  cout << "** Exiting..." << endl;
-  cout<<frac<<endl;
+  globalwatch.Stop();
+
+
+    // Screen output
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**        ################## Time report #################         **"<< endl;
+  cout << left  << setw(32) <<"**              Time report for "<<""
+       << left  << setw(15) << allEntries <<""
+       << right << setw(22) <<"events         **"<<endl;
+  cout <<"**                                                                 **"<< endl;
+  cout << left  << setw(10) <<"**"<<""
+       << left  << setw(15) <<"Time (s):"<<""
+       << right << setw(15) <<"CPU"<<""
+       << right << setw(15) <<"Real"<<""
+       << right << setw(14)  <<"**"<<endl;
+  cout << left  << setw(10) <<"**"<<""
+       << left  << setw(15) <<" +  Global:"<<""
+       << right << setw(15) <<globalwatch.CpuTime()<<""
+       << right << setw(15) <<globalwatch.RealTime()<<""
+       << right << setw(14)  <<"**"<<endl;
+
+
+  string tempstring = detectorcard + " has been used.";
+  cout <<"**        ################## Configuration ###############         **"<< endl;
+  cout << left << setw(16) << "**        " << ""
+       << left << setw(51) <<  tempstring << "**" << endl;
+
+  tempstring =  outputfilename + " has been created.";
+  cout << left << setw(16) << "**        " << ""
+       << left << setw(51) <<  tempstring << "**" << endl;
+
+  cout << left << setw(16) << "**        " << ""
+       << left << setw(51) <<  get_time_date() << "**" << endl;
+
+
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                        Exiting ...                              **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"**                                                                 **"<< endl;
+  cout <<"*********************************************************************"<< endl;
+