Fork me on GitHub

Changeset 467 in svn


Ignore:
Timestamp:
Jul 10, 2009, 5:57:55 PM (15 years ago)
Author:
Xavier Rouby
Message:

cleaning, comments, progressbar, cout

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/CHANGELOG

    r466 r467  
    77   - it is now possible to read "delphes root files" as an input for Delphes (GEN tree)
    88   - bug fix: number of events to process is correct event if multiple files in the inputlist
     9   - updates in Resolution.cpp & Resolution_ATLAS.cpp
    910
    1011 /-----------------------------------------\

  • trunk/Resolutions.cpp

    r458 r467  
    3838#include "TApplication.h"
    3939#include "TFile.h"
     40#include "TStopwatch.h"
    4041
    4142#include "ExRootTreeReader.h"
    4243#include "ExRootTreeWriter.h"
    4344#include "ExRootTreeBranch.h"
     45#include "ExRootProgressBar.h"
    4446#include "TreeClasses.h"
    4547
    … …  
    189191{
    190192  int appargc = 2;
    191   char *appName = "Resolution";
    192   char *appargv[] = {appName, "-b"};
     193  char *appName= new char[20];
     194  char *appOpt= new char[20];
     195  sprintf(appName,"Resolution_ATLAS");
     196  sprintf(appOpt,"-b");
     197  char *appargv[] = {appName,appOpt};
    193198  TApplication app(appName, &appargc, appargv);
    194  
     199  delete [] appName;
     200  delete [] appOpt;
     201
    195202  if(argc != 3) {
    196203      cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
    197       cout << " input_list - list of files in root format," << endl;
     204      cout << " input_file - input file in Delphes-root format," << endl;
    198205      cout << " output_file - output file." << endl;
    199206      exit(1);
    200207  }
    201208
    202   srand (time (NULL));         /* Initialisation du générateur */
    203  
    204   //read the input TROOT file
     209
     210  // 1. ********** initialisation ***********
     211
     212  srand (time (NULL));
     213  TStopwatch globalwatch;
     214  globalwatch.Start();
     215 
    205216  string inputfilename(argv[1]), outputfilename(argv[2]);
    206217
     218  // 1.1 checks the input file
     219  if(inputfilename.find(".root") > inputfilename.length()   ) {
     220    cout << "input_file should be a .root file from Delphes!\n";
     221    return -1;
     222  }
     223  TFile f(inputfilename.c_str());
     224  if (!f.FindKey("GEN") || !f.FindKey("Analysis")           ) {
     225    cout << "input_file should be a .root file from Delphes!\n";
     226    cout << "it should contain both \"GEN\" and \"Analysis\" trees at least.\n";
     227    return -1;
     228  }
     229  f.Close();
     230
     231  // 1.2 checks the output file
    207232  if(outputfilename.find(".root") > outputfilename.length() ) {
    208233    cout << "output_file should be a .root file!\n";
    209234    return -1;
    210235  }
    211 
    212 
    213 
    214236  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
    215237  outputFile->Close();
    216  
     238 
     239  // 1.3 Reads the trees in input file
    217240  TChain chainGEN("GEN");
    218241  chainGEN.Add(inputfilename.c_str());
    … …  
    229252  TIter itGen((TCollection*)branchGen);
    230253 
    231   //write the output root file
     254
     255  // 1.4 Prepares the output root file
    232256  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
    233257  ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
    … …  
    236260  ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
    237261  ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
    238  
     262 
     263  // 1.5 other initialisations
    239264  TRootGenParticle *particle;
    240  
    241265  RESOLELEC * elementElec;
    242266  RESOLMUON *elementMuon;
    … …  
    245269  ETMIS *elementEtmis;
    246270 
    247   int numTau=0;
    248   int numTauRec=0;
     271  int numTau=0, numTauRec=0;
    249272 
    250273  RESOLution *DET = new RESOLution();
     274  /*
     275  string detectorcard = "data/DetectorCard_CMS.dat";
     276  const float dR_jetpairing = 0.25;
     277  const float jet_pt_cut = 1;
     278  */
    251279  DET->ReadDataCard("data/DetectorCard_CMS.dat");
    252280
    253281
    254   //Jet information
    255282  JetsUtil *JETRUN = new JetsUtil("data/DetectorCard_CMS.dat");
    256283
    257   TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
    258   TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information
     284  TLorentzVector genMomentum(0,0,0,0); //generator level information
     285  TLorentzVector recoMomentum(0,0,0,0);//reconstructed level information
    259286  LorentzVector jetMomentum;
    260287 
    261288  vector<fastjet::PseudoJet> input_particlesGEN;//for FastJet algorithm
    262289  vector<fastjet::PseudoJet> sorted_jetsGEN;
    263  
    264290  vector<int> NTrackJet;
    265  
    266291  vector<TLorentzVector> towers;
    267292 
    268   // Loop over all events
     293  // 2. Loop over all events
    269294  Long64_t entry, allEntries = treeReader->GetEntries();
    270   cout << "** Chain contains " << allEntries << " events" << endl;
     295  cout << endl << endl;
     296  cout <<"*********************************************************************"<< endl;
     297  cout <<"**                                                                 **"<< endl;
     298  cout <<"**                                                                 **"<< endl;
     299  cout <<"**                                                                 **"<< endl;
     300  cout <<"**        ####### Start resolution processing    ########          **"<< endl;
     301  cout << left  << setw(52) <<"**              Total number of events to run: "<<""
     302       << left  << setw(15) << allEntries <<""
     303       << right << setw(2) <<"**"<<endl;
     304
     305  ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
     306
    271307  for(entry = 0; entry < allEntries; ++entry)
    272308    {
     309      Progress->Update(entry);
    273310      TLorentzVector PTmisReco(0,0,0,0);
    274311      TLorentzVector PTmisGEN(0,0,0,0);
    … …  
    276313      treeReaderGEN->ReadEntry(entry);
    277314      treeWriter->Clear();
    278       if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
    279315     
    280316      TSimpleArray<TRootGenParticle> bGen;
    … …  
    301337          //Calculate ETMIS from generated particles
    302338          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmisGEN = PTmisGEN + genMomentum;
    303          
     339
     340          //Electrons and muons 
    304341          if( (particle->Status == 1)    &&
    305342            ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
    … …  
    320357                 case pMU: // all muons with eta < DET->MAX_MU
    321358                    PairingMuon(recoMomentum,genMomentum,branchMuon);
    322                     if(recoMomentum.E() !=0){
     359                    if(recoMomentum.E()!=0){
    323360                     elementMuon = (RESOLMUON*) branchmuon->NewEntry();
    324                      elementMuon->OverPT = (1/genMomentum.Pt());
    325                      elementMuon->OverSmearedPT = (1/recoMomentum.Pt());}
     361                     elementMuon->OverPT = 1./genMomentum.Pt();
     362                     elementMuon->OverSmearedPT = 1./recoMomentum.Pt();}
    326363                   break; // case pMU
    327364                 default:
    … …  
    354391      elementEtmis->Ex = (-PTmisGEN).Px();
    355392      elementEtmis->SEt = ScalarEt;
    356 
    357393      elementEtmis->EtSmeare = (PTmisReco).Pt()-(PTmisGEN).Pt();
    358394      elementEtmis->ExSmeare = (-PTmisReco).Px()-(PTmisGEN).Px();
    … …  
    362398
    363399      TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
    364 
    365400      TLorentzVector JETreco(0,0,0,0);
    366401      for (unsigned int i = 0; i < sorted_jetsGEN.size(); i++) {
    … …  
    372407            elementJet= (RESOLJET*) branchjet->NewEntry();
    373408            elementJet->PT = JETgen.Et();
    374             elementJet->SmearedPT = JETreco.Et()/JETgen.Et();
     409            elementJet->SmearedPT = JETreco.Et()/JETgen.Et();  //// la difference pourrait être ici
    375410          }
    376411      }
    377       numTau = numTau+TausHadr.GetEntries();
     412      numTau += TausHadr.GetEntries();
    378413
    379414      TIter itJet((TCollection*)branchJet);
    … …  
    391426                  {
    392427                    elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
    393                     elementTaujet->EnergieCen = (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E());
     428                    elementTaujet->EnergieCen = EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E();
    394429                    elementTaujet->NumTrack = NumTracks(branchTracks,DET->TAU_track_pt,JETT.Eta(),JETT.Phi(),DET->TAU_track_scone);
    395430                    if( (EnergySmallCone(towers,JETT.Eta(),JETT.Phi(),DET->TAU_energy_scone,DET->JET_seed)/JETT.E()) > 0.95
    … …  
    401436       
    402437      } // for itJet : loop on all jets
    403 //cout<<"i"<<endl;
    404438     
    405439      treeWriter->Fill();
    406440    } // Loop over all events
    407441  treeWriter->Write();
    408 float frac = numTauRec/numTau; 
    409 cout<<numTauRec<<endl;
    410 cout<<numTau<<endl;
    411 
    412   cout << "** Exiting..." << endl;
    413   cout<<frac<<endl;
     442  globalwatch.Stop();
     443
     444
     445    // Screen output
     446  cout <<"**                                                                 **"<< endl;
     447  cout <<"**                                                                 **"<< endl;
     448  cout <<"**        ################## Time report #################         **"<< endl;
     449  cout << left  << setw(32) <<"**              Time report for "<<""
     450       << left  << setw(15) << allEntries <<""
     451       << right << setw(22) <<"events         **"<<endl;
     452  cout <<"**                                                                 **"<< endl;
     453  cout << left  << setw(10) <<"**"<<""
     454       << left  << setw(15) <<"Time (s):"<<""
     455       << right << setw(15) <<"CPU"<<""
     456       << right << setw(15) <<"Real"<<""
     457       << right << setw(14)  <<"**"<<endl;
     458  cout << left  << setw(10) <<"**"<<""
     459       << left  << setw(15) <<" +  Global:"<<""
     460       << right << setw(15) <<globalwatch.CpuTime()<<""
     461       << right << setw(15) <<globalwatch.RealTime()<<""
     462       << right << setw(14)  <<"**"<<endl;
     463
     464
     465  string tempstring = detectorcard + " has been used.";
     466  cout <<"**        ################## Configuration ###############         **"<< endl;
     467  cout << left << setw(16) << "**        " << ""
     468       << left << setw(51) <<  tempstring << "**" << endl;
     469
     470  tempstring =  outputfilename + " has been created.";
     471  cout << left << setw(16) << "**        " << ""
     472       << left << setw(51) <<  tempstring << "**" << endl;
     473
     474  cout << left << setw(16) << "**        " << ""
     475       << left << setw(51) <<  get_time_date() << "**" << endl;
     476
     477
     478  cout <<"**                                                                 **"<< endl;
     479  cout <<"**                        Exiting ...                              **"<< endl;
     480  cout <<"**                                                                 **"<< endl;
     481  cout <<"**                                                                 **"<< endl;
     482  cout <<"*********************************************************************"<< endl;
     483
    414484
    415485 

  • trunk/Resolutions_ATLAS.cpp

    r458 r467  
    3838#include "TApplication.h"
    3939#include "TFile.h"
     40#include "TStopwatch.h"
    4041
    4142#include "ExRootTreeReader.h"
    4243#include "ExRootTreeWriter.h"
    4344#include "ExRootTreeBranch.h"
     45#include "ExRootProgressBar.h"
    4446#include "TreeClasses.h"
    4547
    … …  
    189191{
    190192  int appargc = 2;
    191   char *appName = "Resolution";
    192   char *appargv[] = {appName, "-b"};
     193  char *appName= new char[20];
     194  char *appOpt= new char[20];
     195  sprintf(appName,"Resolution_ATLAS");
     196  sprintf(appOpt,"-b");
     197  char *appargv[] = {appName,appOpt};
    193198  TApplication app(appName, &appargc, appargv);
    194  
     199  delete [] appName;
     200  delete [] appOpt;
     201
    195202  if(argc != 3) {
    196203      cout << " Usage: " << argv[0] << " input_file" << " output_file" << endl;
    … …  
    200207  }
    201208
    202   srand (time (NULL));         /* Initialisation du générateur */
     209
     210  // 1. ********** initialisation ***********
     211
     212  srand (time (NULL));         
     213  TStopwatch globalwatch;
     214  globalwatch.Start();
    203215 
    204   //read the input TROOT file
    205216  string inputfilename(argv[1]), outputfilename(argv[2]);
    206217
     218  // 1.1 checks the input file
     219  if(inputfilename.find(".root") > inputfilename.length()   ) {
     220    cout << "input_file should be a .root file from Delphes!\n";
     221    return -1;
     222  }
     223  TFile f(inputfilename.c_str());
     224  if (!f.FindKey("GEN") || !f.FindKey("Analysis")           ) {
     225    cout << "input_file should be a .root file from Delphes!\n";
     226    cout << "it should contain both \"GEN\" and \"Analysis\" trees at least.\n";
     227    return -1;
     228  }
     229  f.Close();
     230
     231  // 1.2 checks the output file
    207232  if(outputfilename.find(".root") > outputfilename.length() ) {
    208233    cout << "output_file should be a .root file!\n";
    209234    return -1;
    210235  }
    211 
    212 
    213 
    214236  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE");// Creates the file, but should be closed just after
    215237  outputFile->Close();
    216  
     238 
     239  // 1.3 Reads the trees in  input file
    217240  TChain chainGEN("GEN");
    218241  chainGEN.Add(inputfilename.c_str());
    … …  
    228251  const TClonesArray *branchGen = treeReaderGEN->UseBranch("Particle");
    229252  TIter itGen((TCollection*)branchGen);
    230  
    231   //write the output root file
     253 
     254 
     255  // 1.4 Prepares the output root file
    232256  ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputfilename, "Analysis");
    233257  ExRootTreeBranch *branchjet = treeWriter->NewBranch("JetPTResol", RESOLJET::Class());
    … …  
    236260  ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
    237261  ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
    238  
     262 
     263  // 1.5 other initialisations
    239264  TRootGenParticle *particle;
    240  
    241265  RESOLELEC * elementElec;
    242266  RESOLMUON *elementMuon;
    … …  
    245269  ETMIS *elementEtmis;
    246270 
    247   int numTau=0;
    248   int numTauRec=0;
     271  int numTau=0, numTauRec=0;
    249272 
    250273  RESOLution *DET = new RESOLution();
    … …  
    259282  DET->ReadDataCard(detectorcard);
    260283
    261 
    262   //Jet information
    263284  JetsUtil *JETRUN = new JetsUtil(detectorcard);
    264285
    265   TLorentzVector genMomentum(0,0,0,0);//TLorentzVector containing generator level information
    266   TLorentzVector recoMomentum(0,0,0,0);//TLorentzVector containing generator level information
     286  TLorentzVector genMomentum(0,0,0,0); // generator level information
     287  TLorentzVector recoMomentum(0,0,0,0);// reconstruction level information
    267288  LorentzVector jetMomentum;
    268289 
    … …  
    272293  vector<TLorentzVector> towers;
    273294 
    274   // Loop over all events
     295  // 2. Loop over all events
    275296  Long64_t entry, allEntries = treeReader->GetEntries();
    276   cout << "** Chain contains " << allEntries << " events" << endl;
     297  cout << endl << endl;
     298  cout <<"*********************************************************************"<< endl;
     299  cout <<"**                                                                 **"<< endl;
     300  cout <<"**                                                                 **"<< endl;
     301  cout <<"**                                                                 **"<< endl;
     302  cout <<"**        ####### Start resolution processing    ########          **"<< endl;
     303  cout << left  << setw(52) <<"**              Total number of events to run: "<<""
     304       << left  << setw(15) << allEntries <<""
     305       << right << setw(2) <<"**"<<endl;
     306
     307  ExRootProgressBar *Progress = new ExRootProgressBar(allEntries);
     308
    277309  for(entry = 0; entry < allEntries; ++entry)
    278310    {
     311      Progress->Update(entry);
    279312      TLorentzVector PTmisReco(0,0,0,0);
    280313      TLorentzVector PTmisGEN(0,0,0,0);
    … …  
    282315      treeReaderGEN->ReadEntry(entry);
    283316      treeWriter->Clear();
    284       if((entry % 100) == 0 && entry > 0 )  cout << "** Processing element # " << entry << endl;
    285317     
    286318      TSimpleArray<TRootGenParticle> bGen;
    … …  
    372404        TLorentzVector JETgen(0,0,0,0);
    373405        JETgen.SetPxPyPzE(sorted_jetsGEN[i].px(),sorted_jetsGEN[i].py(),sorted_jetsGEN[i].pz(),sorted_jetsGEN[i].E());
    374         PairingJet(JETreco,JETgen,branchJet);
     406        PairingJet(JETreco,JETgen,branchJet,dR_jetpairing);
    375407        if(JETreco.Pt()>jet_pt_cut)
    376408          {
    … …  
    383415          }
    384416      }
    385       numTau = numTau+TausHadr.GetEntries();
     417      numTau += TausHadr.GetEntries();
    386418
    387419      TIter itJet((TCollection*)branchJet);
    … …  
    413445    } // Loop over all events
    414446  treeWriter->Write();
    415 
    416   cout << detectorcard << " has been used.\n";
    417   cout << "** Exiting..." << endl;
     447  globalwatch.Stop();
     448
     449
     450  // Screen output
     451  cout <<"**                                                                 **"<< endl;
     452  cout <<"**                                                                 **"<< endl;
     453  cout <<"**        ################## Time report #################         **"<< endl;
     454  cout << left  << setw(32) <<"**              Time report for "<<""
     455       << left  << setw(15) << allEntries <<""
     456       << right << setw(22) <<"events         **"<<endl;
     457  cout <<"**                                                                 **"<< endl;
     458  cout << left  << setw(10) <<"**"<<""
     459       << left  << setw(15) <<"Time (s):"<<""
     460       << right << setw(15) <<"CPU"<<""
     461       << right << setw(15) <<"Real"<<""
     462       << right << setw(14)  <<"**"<<endl;
     463  cout << left  << setw(10) <<"**"<<""
     464       << left  << setw(15) <<" +  Global:"<<""
     465       << right << setw(15) <<globalwatch.CpuTime()<<""
     466       << right << setw(15) <<globalwatch.RealTime()<<""
     467       << right << setw(14)  <<"**"<<endl;
     468
     469
     470  string tempstring = detectorcard + " has been used.";
     471  cout <<"**        ################## Configuration ###############         **"<< endl;
     472  cout << left << setw(16) << "**        " << ""
     473       << left << setw(51) <<  tempstring << "**" << endl;
     474
     475  tempstring =  outputfilename + " has been created.";
     476  cout << left << setw(16) << "**        " << ""
     477       << left << setw(51) <<  tempstring << "**" << endl;
     478
     479  cout << left << setw(16) << "**        " << ""
     480       << left << setw(51) <<  get_time_date() << "**" << endl;
     481
     482  cout <<"**                                                                 **"<< endl;
     483  cout <<"**                        Exiting ...                              **"<< endl;
     484  cout <<"**                                                                 **"<< endl;
     485  cout <<"**                                                                 **"<< endl;
     486  cout <<"*********************************************************************"<< endl;
     487
    418488 
    419489  delete treeWriter;
    … …  
    422492}
    423493
     494
     495
     496
     497
     498
     499
     500
     501
     502
     503
     504
     505
     506
     507
     508
     509
     510
     511
     512
     513
     514
     515
     516
Note: See TracChangeset for help on using the changeset viewer.