Changeset 94 in svn for trunk/Delphes.cpp

Timestamp:

Dec 12, 2008, 5:32:29 PM (16 years ago)

Author:

severine ovyn

Message:

Add frog plus cleaning + bugs remove

File:

• trunk/Delphes.cpp (modified) (12 diffs)

trunk/Delphes.cpp

@@ -29 +29 @@
 #include "interface/VeryForward.h"
 #include "interface/JetUtils.h"
+#include "interface/FrogUtil.h"
 
 #include <vector>
@@ -40 +41 @@
   cout << "** TODO list ..." << endl;
   while(infile.good()) {
-        string temp;
-        getline(infile,temp);
-        cout << "*" << temp << endl;
+    string temp;
+    getline(infile,temp);
+    cout << "*" << temp << endl;
   }
   cout << "** done...\n"; 
@@ -55 +56 @@
   char *appargv[] = {appName, "-b"};
   TApplication app(appName, &appargc, appargv);
-
+  
   if(argc != 4 && argc != 3 && argc != 5) {
-      cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl;
-      cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
-      cout << " output_file - output file." << endl;
-      cout << " detector_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
-      cout << " trigger_card - Datacard containing the trigger algorithms (optional) "<<endl;
-      exit(1);
+    cout << " Usage: " << argv[0] << " input_file output_file [detector_card] [trigger_card] " << endl;
+    cout << " input_list - list of files in Ntpl, StdHep of LHEF format," << endl;
+    cout << " output_file - output file." << endl;
+    cout << " detector_card - Datacard containing resolution variables for the detector simulation (optional) "<<endl;
+    cout << " trigger_card - Datacard containing the trigger algorithms (optional) "<<endl;
+    exit(1);
   }
-
+  
   srand (time (NULL));         /* Initialisation du générateur */
   
@@ -70 +71 @@
   string inputFileList(argv[1]), outputfilename(argv[2]);
   if(outputfilename.find(".root") > outputfilename.length() ) {
-        cout << "output_file should be a .root file!\n"; 
-        exit(1);
+    cout << "output_file should be a .root file!\n"; 
+    exit(1);
   }
   //create output log-file name
@@ -77 +78 @@
   string LogName = forLog.erase(forLog.find(".root"));
   LogName = LogName+"_run.log";
-
+  
   TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
   outputFile->Close();
-
+  
   string line;
   ifstream infile(inputFileList.c_str());
   infile >> line; // the first line determines the type of input files
-
+  
   //read the datacard input file
   string DetDatacard("");
-  if(argc==4)  DetDatacard =argv[3];
-
+  if(argc>=4)  DetDatacard =argv[3];
+  
   //Smearing information
   RESOLution *DET = new RESOLution();
   DET->ReadDataCard(DetDatacard);
   DET->Logfile(LogName);
-
+  
   //read the trigger input file
   string TrigDatacard("data/trigger.dat");
   if(argc==5)  TrigDatacard =argv[4];
- 
+  
   //Trigger information
   TriggerTable *TRIGT = new TriggerTable();
   TRIGT->TriggerCardReader(TrigDatacard.c_str());
   TRIGT->PrintTriggerTable(LogName);
-
+  
   //Propagation of tracks in the B field
   TrackPropagation *TRACP = new TrackPropagation(); 
-
+  
   //Jet information
   JetsUtil *JETRUN = new JetsUtil();
-
+  
   //VFD information
   VeryForward * VFD = new VeryForward();
-
+  
   //todo(LogName.c_str());
- 
+  
   DataConverter *converter=0;
-   
+  
   if(strstr(line.c_str(),".hep"))
     {                           
@@ -174 +175 @@
   TLorentzVector genMomentumCalo(0,0,0,0);
   LorentzVector jetMomentum;
- 
+  
   vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
   vector<fastjet::PseudoJet> sorted_jets;
-
+  
   vector<TLorentzVector> TrackCentral;  
   vector<PhysicsTower> towers;
@@ -186 +187 @@
   TSimpleArray<TRootGenParticle> NFCentralQ;
   
-
-
+  
+  
   // Loop over all events
   Long64_t entry, allEntries = treeReader->GetEntries();
@@ -212 +213 @@
         {
           int pid = abs(particle->PID);
-         //// This subarray is needed for the B-jet algorithm
+          //// This subarray is needed for the B-jet algorithm
           // optimization for speed : put first PID condition, then ETA condition, then either pt or status
           if( (pid <= pB || pid == pGLUON) &&// is it a light quark or a gluon, i.e. is it one of these : u,d,c,s,b,g ?
-              fabs(particle->Eta) < DET->MAX_TRACKER &&
+              fabs(particle->Eta) < DET->CEN_max_tracker &&
               particle->Status != 1 &&
               particle->PT > DET->PT_QUARKS_MIN ) {
             NFCentralQ.Add(particle);
           }
-                  
+          
           // keeps only final particles, visible by the central detector, including the fiducial volume
           // the ordering of conditions have been optimised for speed : put first the STATUS condition
@@ -227 +228 @@
           if( (particle->Status == 1)    &&
               ((pid != pNU1) && (pid != pNU2) && (pid != pNU3)) &&
-              (fabs(particle->Eta) < DET->MAX_CALO_FWD)
-            ) 
+              (fabs(particle->Eta) < DET->CEN_max_calo_fwd)
+              ) 
             { 
-            genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
-            TRACP->Propagation(particle,genMomentum);
-            float eta=fabs(genMomentum.Eta());
-
-            switch(pid) {
+              genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
+              if(DET->FLAG_bfield==1)TRACP->Propagation(particle,genMomentum);
+              float eta=fabs(genMomentum.Eta());
               
-            case pE: // all electrons with eta < DET->MAX_CALO_FWD
-              DET->SmearElectron(genMomentum);
-              if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER && genMomentum.Pt() > DET->ELEC_pt){
-              electron.push_back(ParticleUtil(genMomentum,particle->PID));
-              }
-              break; // case pE
-            case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
-              DET->SmearElectron(genMomentum);
-              if(genMomentum.E()!=0 && eta < DET->MAX_TRACKER && genMomentum.Pt() > DET->GAMMA_pt) {
-                gamma.push_back(ParticleUtil(genMomentum,particle->PID));
+              switch(pid) {
+                
+              case pE: // all electrons with eta < DET->MAX_CALO_FWD
+                DET->SmearElectron(genMomentum);
+                if(genMomentum.E()!=0 && eta < DET->CEN_max_tracker && genMomentum.Pt() > DET->PTCUT_elec){
+                  electron.push_back(ParticleUtil(genMomentum,particle->PID));
+                }
+                break; // case pE
+              case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
+                DET->SmearElectron(genMomentum);
+                if(genMomentum.E()!=0 && eta < DET->CEN_max_tracker && genMomentum.Pt() > DET->PTCUT_gamma) {
+                  gamma.push_back(ParticleUtil(genMomentum,particle->PID));
+                }
+                break; // case pGAMMA
+              case pMU: // all muons with eta < DET->MAX_MU
+                DET->SmearMu(genMomentum);
+                if(genMomentum.E()!=0 && eta < DET->CEN_max_mu &&  genMomentum.Pt() > DET->PTCUT_muon){
+                  muon.push_back(ParticleUtil(genMomentum,particle->PID));
+                }
+                break; // case pMU
+              case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
+              case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
+                DET->SmearHadron(genMomentum, 0.7);
+                break; // case hadron
+              default:   // all other final particles with eta < DET->MAX_CALO_FWD
+                DET->SmearHadron(genMomentum, 1.0);
+                break;
+              } // switch (pid)
+              
+              // all final particles but muons and neutrinos    
+              // for calorimetric towers and mission PT
+              int charge=Charge(pid); 
+              if(genMomentum.E() !=0 && pid != pMU) {
+                if(charge == 0 || (charge !=0 && genMomentum.Pt() >= DET->TRACK_ptmin)){
+                  PhysicsTower CaloTower = PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+                  towers.push_back(CaloTower);
+                  // create a fastjet::PseudoJet with these components and put it onto
+                  // back of the input_particles vector
+                  input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
+                  
+                  genMomentumCalo.SetPxPyPzE(CaloTower.fourVector.px,CaloTower.fourVector.py,CaloTower.fourVector.pz,CaloTower.fourVector.E);
+                  
+                  elementCalo = (TRootCalo*) branchCalo->NewEntry();
+                  elementCalo->Set(genMomentumCalo);
+                  DET->BinEtaPhi(genMomentumCalo.Phi(), genMomentumCalo.Eta(), elementCalo->Phi, elementCalo->Eta);
+                }
               }
-              break; // case pGAMMA
-            case pMU: // all muons with eta < DET->MAX_MU
-              DET->SmearMu(genMomentum);
-              if(genMomentum.E()!=0 && eta < DET->MAX_MU &&  genMomentum.Pt() > DET->MUON_pt){
-              muon.push_back(ParticleUtil(genMomentum,particle->PID));
-              }
-              break; // case pMU
-            case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
-            case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
-              DET->SmearHadron(genMomentum, 0.7);
-              break; // case hadron
-            default:   // all other final particles with eta < DET->MAX_CALO_FWD
-              DET->SmearHadron(genMomentum, 1.0);
-              break;
-            } // switch (pid)
-            
-            // all final particles but muons and neutrinos      
-            // for calorimetric towers and mission PT
-            int charge=Charge(pid); 
-            if(genMomentum.E() !=0 && pid != pMU) {
-              if(charge == 0 || (charge !=0 && genMomentum.Pt() >= DET->PT_TRACKS_MIN)){
-                PhysicsTower CaloTower = PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
-                towers.push_back(CaloTower);
-                // create a fastjet::PseudoJet with these components and put it onto
-                // back of the input_particles vector
-                input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())); 
-                
-                genMomentumCalo.SetPxPyPzE(CaloTower.fourVector.px,CaloTower.fourVector.py,CaloTower.fourVector.pz,CaloTower.fourVector.E);
-                
-                elementCalo = (TRootCalo*) branchCalo->NewEntry();
-                elementCalo->Set(genMomentumCalo);
-                DET->BinEtaPhi(genMomentumCalo.Phi(), genMomentumCalo.Eta(), elementCalo->Phi, elementCalo->Eta);
-              }
-            }
-            
-            // all final charged particles
-            if(
-               (genMomentum.E()!=0) &&
-               (fabs(genMomentum.Eta()) < DET->MAX_TRACKER) && 
-               (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
-               ((rand()%100) < DET->TRACKING_EFF)  &&
-               (charge!=0)
-               )
-              {
-                elementTracks = (TRootTracks*) branchTracks->NewEntry();
-                elementTracks->Set(genMomentum);
-                TrackCentral.push_back(genMomentum);
-              }
-            
+              
+              // all final charged particles
+              if(
+                 (genMomentum.E()!=0) &&
+                 (fabs(genMomentum.Eta()) < DET->CEN_max_tracker) && 
+                 (genMomentum.Pt() > DET->TRACK_ptmin ) &&     // pt too small to be taken into account
+                 ((rand()%100) < DET->TRACK_eff)  &&
+                 (charge!=0)
+                 )
+                {
+                  elementTracks = (TRootTracks*) branchTracks->NewEntry();
+                  elementTracks->Set(genMomentum);
+                  TrackCentral.push_back(genMomentum);
+                }
+              
             } // switch
           
-          VFD->ZDC(treeWriter,branchZDC,particle);
-          VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
+          if(DET->FLAG_vfd==1)
+            {  
+              VFD->ZDC(treeWriter,branchZDC,particle);
+              VFD->RomanPots(treeWriter,branchRP220,branchFP420,particle);
+            }
           
         } // while
@@ -350 +354 @@
   treeWriter->Write();
   delete treeWriter;
-
+  
   //running the trigger in case the FLAG trigger is put to 1 in the datacard
- 
-  if(DET->DOTRIGGER == 1)
+  
+  if(DET->FLAG_trigger == 1)
     {
       TChain chainT("Analysis");
@@ -378 +382 @@
           treeWriterT->Fill();
         }
-
+      
       treeWriterT->Write();
       delete treeWriterT;
     }
-
+  
+  //FROG display
+  if(DET->FLAG_frog == 1)
+    {
+      FrogDisplay *FROG = new FrogDisplay();
+      FROG->BuidEvents(outputfilename,DET->NEvents_Frog);
+      FROG->BuildGeom();
+    }
+  
   cout << "** Exiting..." << endl;
   
@@ -391 +403 @@
   delete JETRUN;
   delete VFD;
- 
+  
   if(converter) delete converter;
-
+  
   todo("TODO");
 }