Fork me on GitHub

Changeset 26 in svn


Ignore:
Timestamp:
Nov 13, 2008, 9:02:15 AM (16 years ago)
Author:
severine ovyn
Message:

OK for jet resolution

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Resolutions.cpp

    r24 r26  
    5252//------------------------------------------------------------------------------
    5353
    54 void PairingJet(TLorentzVector &JETSm, const TLorentzVector& JET,  vector<fastjet::PseudoJet> sorted_jetsS)
     54
     55
     56void PairingJet(TLorentzVector &JETSm, TLorentzVector JET,  vector<fastjet::PseudoJet> sorted_jetsS)
    5557{
    5658  JETSm.SetPxPyPzE(0,0,0,0);
    … …  
    7072}
    7173
    72 //------------------------------------------------------------------------------
    7374
    7475int main(int argc, char *argv[])
    7576{
    7677  int appargc = 2;
    77   char appName[100]; sprintf(appName,argv[0]);
     78  char *appName = "Smearing";
    7879  char *appargv[] = {appName, "-b"};
    7980  TApplication app(appName, &appargc, appargv);
    … …  
    9798  TFile *outputFile = TFile::Open(outputfilename.c_str(), "RECREATE"); // Creates the file, but should be closed just after
    9899  outputFile->Close();
    99 
     100 
    100101  string line;
    101102  ifstream infile(inputFileList.c_str());
    102103  infile >> line; // the first line determines the type of input files
    103  
     104  
    104105  DataConverter *converter=0;
    105    
     106 
    106107  if(strstr(line.c_str(),".hep"))
    107108    {                           
    … …  
    129130    }
    130131  else { cout << "***  " << line.c_str() << "\n***  file format not identified\n***  Exiting\n"; return -1;};
    131 
     132 
    132133  TChain chain("GEN");
    133134  chain.Add(outputfilename.c_str());
    … …  
    143144  ExRootTreeBranch *branchtaujet = treeWriter->NewBranch("TauJetPTResol", TAUHAD::Class());
    144145  ExRootTreeBranch *branchetmis = treeWriter->NewBranch("ETmisResol",ETMIS::Class());
    145 
     146 
    146147  TRootGenParticle *particle;
    147148  TRootETmis *etmisc;
    148 
     149 
    149150  RESOLELEC *elementElec;
    150151  RESOLMUON *elementMuon;
    … …  
    152153  TAUHAD *elementTaujet;
    153154  ETMIS *elementEtmis;
    154 
    155  
     155 
     156  
    156157  //read the datacard input file
    157158  string DetDatacard("");
    … …  
    163164  LorentzVector jetMomentum;
    164165  vector<TLorentzVector> TrackCentral;
    165 
     166 
    166167  vector<fastjet::PseudoJet> input_particles;//for FastJet algorithm
    167168  vector<fastjet::PseudoJet> inclusive_jets;
    168169  vector<fastjet::PseudoJet> sorted_jets;
    169 
     170 
    170171  vector<fastjet::PseudoJet> input_particlesS;//for FastJet algorithm
    171172  vector<fastjet::PseudoJet> inclusive_jetsS;
    172173  vector<fastjet::PseudoJet> sorted_jetsS;
    173174  vector<PhysicsTower> towers;
    174 
     175 
    175176  fastjet::JetDefinition jet_def;
    176177  fastjet::JetDefinition jet_defS;
    177178  fastjet::JetDefinition::Plugin * plugins;
    178179  fastjet::JetDefinition::Plugin * pluginsS;
    179 
     180 
    180181  // set up a CDF midpoint jet definition
    181   #ifdef ENABLE_PLUGIN_CDFCONES
     182#ifdef ENABLE_PLUGIN_CDFCONES
    182183  plugins = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    183184  jet_def = fastjet::JetDefinition(plugins);
    184   #else
     185#else
    185186  plugins = NULL;
    186   #endif
    187 
     187#endif
     188 
    188189  // set up a CDF midpoint jet definition
    189   #ifdef ENABLE_PLUGIN_CDFCONES
    190   pluginsS = new fastjet::CDFJetCluPlugin(2,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
     190#ifdef ENABLE_PLUGIN_CDFCONES
     191  pluginsS = new fastjet::CDFJetCluPlugin(1,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
    191192  jet_defS = fastjet::JetDefinition(pluginsS);
    192   #else
     193#else
    193194  pluginsS = NULL;
    194   #endif
    195 
    196  
     195#endif
     196 
     197  
    197198  // Loop over all events
    198199  Long64_t entry, allEntries = treeReader->GetEntries();
    … …  
    218219      sorted_jetsS.clear();
    219220      towers.clear();
    220 
     221     
    221222      // Loop over all particles in event
    222223      while( (particle = (TRootGenParticle*) itGen.Next()) )
    223224        {
    224225          genMomentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
    225 
     226         
    226227          int pid = abs(particle->PID);
    227228          float eta = fabs(particle->Eta);
    228 
     229         
    229230          if(particle->Status == 1)
    230              {
    231                input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    232              }
    233 
     231            {
     232              input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     233            }
     234         
    234235          // keeps only final particles, visible by the central detector, including the fiducial volume
    235236          // the ordering of conditions have been optimised for speed : put first the STATUS condition
    236237          if( (particle->Status == 1)    &&
    237238              (
    238                 (pid == pMU  && eta < DET->MAX_MU) ||
    239                 (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD)
    240               )
    241             ) {
    242                 if(pid != pMU)PTmis = PTmis + genMomentum;//ptmis
    243                 switch(pid) {
    244                    
    245                     case pE: // all electrons with eta < DET->MAX_CALO_FWD
    246                       DET->SmearElectron(genMomentum);
    247                       break; // case pE
    248 
    249                     case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
    250                       DET->SmearElectron(genMomentum);
    251                       break; // case pGAMMA
    252 
    253                     case pMU: // all muons with eta < DET->MAX_MU
    254                       DET->SmearMu(genMomentum);
    255                       break; // case pMU
    256 
    257                     case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
    258                     case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
    259                       DET->SmearHadron(genMomentum, 0.7);
    260                       break; // case hadron
    261                      
    262                     default:   // all other final particles with eta < DET->MAX_CALO_FWD
    263                       DET->SmearHadron(genMomentum, 1.0);
    264                       break;
    265                 } // switch (pid)
    266 
    267                 // all final particles but muons and neutrinos 
    268                 // for calorimetric towers and mission PT
    269                 //if(genMomentum.E()!=0) PTmis = PTmis + genMomentum;//ptmis
    270 
    271                 if(pid != pMU)
    272                   {
    273                      PTmisS = PTmisS + genMomentum;
    274                      towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
    275                      input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
    276                   }
    277                
    278                 // all final charged particles
    279                 if(
    280                      ((rand()%100) < DET->TRACKING_EFF)  &&
    281                      (genMomentum.E()!=0) &&
    282                      (fabs(particle->Eta) < DET->MAX_TRACKER) &&
    283                      (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
    284                      (pid != pGAMMA)    &&
    285                      (pid != pPI0)      &&
    286                      (pid != pK0L)      &&
    287                      (pid != pN)        &&
    288                      (pid != pSIGMA0)   &&
    289                      (pid != pDELTA0)   &&
    290                      (pid != pK0S)   // not charged particles : invisible by tracker
     239               (pid == pMU  && eta < DET->MAX_MU) ||
     240               (pid != pMU  && (pid != pNU1) && (pid != pNU2) && (pid != pNU3) && eta < DET->MAX_CALO_FWD) )
     241              )
     242            {
     243              PTmis = PTmis + genMomentum;//ptmis
     244              switch(pid) {
     245
     246              case pE: // all electrons with eta < DET->MAX_CALO_FWD
     247                DET->SmearElectron(genMomentum);
     248                break; // case pE
     249              case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
     250                DET->SmearElectron(genMomentum);
     251                break; // case pGAMMA
     252              case pMU: // all muons with eta < DET->MAX_MU
     253                DET->SmearMu(genMomentum);
     254                break; // case pMU
     255              case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
     256              case pK0S:    // all K0s with eta < DET->MAX_CALO_FWD
     257                DET->SmearHadron(genMomentum, 0.7);
     258                break; // case hadron
     259              default:   // all other final particles with eta < DET->MAX_CALO_FWD
     260                DET->SmearHadron(genMomentum, 1.0);
     261                break;
     262              } // switch (pid)
     263             
     264              if(pid != pMU)
     265                {
     266//                PTmisS = PTmisS + genMomentum;
     267                  towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
     268                  if(genMomentum.Et()>0.5)input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
     269                }
     270             
     271              // all final charged particles
     272              if(
     273                 ((rand()%100) < DET->TRACKING_EFF)  &&
     274                 (genMomentum.E()!=0) &&
     275                 (fabs(particle->Eta) < DET->MAX_TRACKER) &&
     276                 (genMomentum.Pt() > DET->PT_TRACKS_MIN ) &&     // pt too small to be taken into account
     277                 (pid != pGAMMA)    &&
     278                 (pid != pPI0)      &&
     279                 (pid != pK0L)      &&
     280                 (pid != pN)        &&
     281                 (pid != pSIGMA0)   &&
     282                 (pid != pDELTA0)   &&
     283                 (pid != pK0S)   // not charged particles : invisible by tracker
    291284                 )
    292                       TrackCentral.push_back(genMomentum);
     285                TrackCentral.push_back(genMomentum);
    293286            } // switch
    294287        } // while
    295 
    296   TLorentzVector toWerS(0,0,0,0);
    297   //calcul de ETmis au depart des tours calorimetriques..
    298 /*  for(unsigned int i=0; i < towers.size(); i++) {
    299       if(towers[i].fourVector.pt() < 0.5) continue;
    300           toWerS.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
    301           PTmisS = PTmisS + toWerS;
    302           }
    303 */
    304 //PTmis = (0,0,0,0)-PTmis;
    305 //PTmisS = (0,0,0,0)-PTmisS;
    306 
    307 //float ET=PTmis.Et()<<endl;
    308 //float ETS=PTmisS.Et()<<endl;
    309 
    310 cout<<"Ptmis "<<PTmis.Et()<<" PTmis smeare "<<PTmisS.Et()<<endl;
    311 cout<<"Ptmis "<<(-PTmis).Pt()<<" PTmis smeare "<<(-PTmisS).Pt()<<endl;
    312 
    313   elementEtmis= (ETMIS*) branchetmis->NewEntry();
    314   elementEtmis->Et = (PTmis).Et();
    315   elementEtmis->EtSmeare = (PTmisS).Et();
    316 
    317 
    318 
     288     
     289     TLorentzVector TowerTLV(0.,0.,0.,0.);
     290     for(unsigned int i=0; i < towers.size(); i++)
     291       {
     292         TowerTLV.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
     293         if(TowerTLV.Et()>0.5)PTmisS = PTmisS + TowerTLV;
     294      }
     295
     296 
     297      elementEtmis= (ETMIS*) branchetmis->NewEntry();
     298      elementEtmis->Et = (PTmis).Et();
     299      elementEtmis->EtSmeare = (PTmisS).Et()-(PTmis).Et();
     300     
    319301      //*****************************
    320 
    321      double ptmin=1;
    322      if(input_particles.size()!=0)
     302     
     303      double ptmin=1;
     304      if(input_particles.size()!=0)
    323305        {
    324306          fastjet::ClusterSequence clust_seq(input_particles, jet_def);
    … …  
    326308          sorted_jets = sorted_by_pt(inclusive_jets);
    327309        }
    328 
    329      if(input_particlesS.size()!=0)
     310     
     311      if(input_particlesS.size()!=0)
    330312        {
    331313          fastjet::ClusterSequence clust_seqS(input_particlesS, jet_defS);
    … …  
    333315          sorted_jetsS = sorted_by_pt(inclusive_jetsS);
    334316        }
    335 
    336     TLorentzVector JETSm(0,0,0,0);
    337     for (unsigned int i = 0; i < sorted_jets.size(); i++) {
     317     
     318      TLorentzVector JETSm(0,0,0,0);
     319      for (unsigned int i = 0; i < sorted_jets.size(); i++) {
    338320        TLorentzVector JET(0,0,0,0);
    339321        JET.SetPxPyPzE(sorted_jets[i].px(),sorted_jets[i].py(),sorted_jets[i].pz(),sorted_jets[i].E());
    … …  
    341323        if(JETSm.Pt()>1)
    342324          {
    343             float NonSmeareEt=(JET.E()*sin(JET.Theta()));
    344             float SmeareEt=(JETSm.E()*sin(JETSm.Theta()));
    345 
    346325            elementJet= (RESOLJET*) branchjet->NewEntry();
    347326            elementJet->NonSmearePT = JET.Et();
    348327            elementJet->SmearePT = JETSm.Et()/JET.Et();
    349          
     328           
    350329          }
    351 
    352         } // for itJet : loop on all jets
    353 
     330       
     331      } // for itJet : loop on all jets
     332     
    354333      treeWriter->Fill();
    355334    } // Loop over all events
Note: See TracChangeset for help on using the changeset viewer.