Changeset 39 in svn for trunk/Resolutions.cpp

Timestamp:

Nov 18, 2008, 10:30:58 AM (16 years ago)

Author:

severine ovyn

Message:

ok for ETmis

File:

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

trunk/Resolutions.cpp

@@ -52 +52 @@
 //------------------------------------------------------------------------------
 
-
-
-void PairingJet(TLorentzVector &JETSm, TLorentzVector JET,  vector<fastjet::PseudoJet> sorted_jetsS)
+// //********************************** PYTHIA INFORMATION*********************************
+
+TSimpleArray<TRootGenParticle> TauHadr(const TClonesArray *GEN)
+ {
+   TIter it((TCollection*)GEN);
+   it.Reset();
+   TRootGenParticle *gen1;
+   TSimpleArray<TRootGenParticle> array,array2;
+
+   while((gen1 = (TRootGenParticle*) it.Next()))
+      {
+        array.Add(gen1);
+      }
+   it.Reset();
+   bool tauhad;
+   while((gen1 = (TRootGenParticle*) it.Next()))
+    {
+       tauhad=false;
+       if(abs(gen1->PID)==15)
+         {
+            int d1=gen1->D1;
+            int d2=gen1->D2;
+
+            if((d1 < array.GetEntries()) && (d1 > 0) && (d2 < array.GetEntries()) && (d2 > 0))
+               {
+                 tauhad=true;
+                 for(int d=d1; d < d2+1; d++)
+                    {
+                      if(abs(array[d]->PID)== pE || abs(array[d]->PID)== pMU)tauhad=false;
+                    }
+               }
+         }
+        if(tauhad)array2.Add(gen1);
+     }
+   return array2;
+ }
+
+
+
+void PairingJet(TLorentzVector &JETSm, const TLorentzVector &JET,  vector<fastjet::PseudoJet> sorted_jetsS)
 {
   JETSm.SetPxPyPzE(0,0,0,0);
@@ -146 +183 @@
   
   TRootGenParticle *particle;
-  TRootETmis *etmisc;
-  
-  RESOLELEC *elementElec;
+  
+  RESOLELEC * elementElec;
   RESOLMUON *elementMuon;
   RESOLJET *elementJet;
@@ -172 +208 @@
   vector<fastjet::PseudoJet> inclusive_jetsS;
   vector<fastjet::PseudoJet> sorted_jetsS;
+
+  vector<fastjet::PseudoJet> input_particlesT;//for FastJet algorithm
+  vector<fastjet::PseudoJet> inclusive_jetsT;
+  vector<fastjet::PseudoJet> sorted_jetsT;
+
   vector<PhysicsTower> towers;
   
   fastjet::JetDefinition jet_def;
   fastjet::JetDefinition jet_defS;
+  fastjet::JetDefinition jet_defT;
   fastjet::JetDefinition::Plugin * plugins;
   fastjet::JetDefinition::Plugin * pluginsS;
+  fastjet::JetDefinition::Plugin * pluginsT;
   
   // set up a CDF midpoint jet definition
@@ -195 +238 @@
 #endif
   
+ // set up a CDF midpoint jet definition
+ #ifdef ENABLE_PLUGIN_CDFCONES
+ pluginsT = new fastjet::CDFJetCluPlugin(0,DET->CONERADIUS,DET->C_ADJACENCYCUT,DET->C_MAXITERATIONS,DET->C_IRATCH,DET->C_OVERLAPTHRESHOLD);
+ jet_defT = fastjet::JetDefinition(pluginsT);
+ #else
+ pluginsT = NULL;
+ #endif
+ 
   
   // Loop over all events
@@ -212 +263 @@
       TrackCentral.clear();
       TSimpleArray<TRootGenParticle> NFCentralQ;
+   
+      sorted_jetsS.clear();
+      input_particlesS.clear();
+      inclusive_jetsS.clear();
+
+      sorted_jetsT.clear();
+      input_particlesT.clear();
+      inclusive_jetsT.clear();
+
+      sorted_jets.clear();
       input_particles.clear();
       inclusive_jets.clear();
-      sorted_jets.clear();
-      input_particlesS.clear();
-      inclusive_jetsS.clear();
-      sorted_jetsS.clear();
       towers.clear();
       
@@ -228 +285 @@
           float eta = fabs(particle->Eta);
           
-          if(particle->Status == 1)
+          if(particle->Status == 1 && eta < DET->MAX_CALO_FWD)
             {
               input_particles.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
             }
           
+          if((pid == pNU1) || (pid == pNU2) || (pid == pNU3))PTmis = PTmis + genMomentum;
           // 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
@@ -241 +299 @@
               ) 
             {
-              PTmis = PTmis + genMomentum;//ptmis
+              //if((pid != pNU1) && (pid != pNU2) && (pid != pNU3))PTmis = PTmis + genMomentum;//ptmis
+
+              Float_t nonS=0;
               switch(pid) {
-
               case pE: // all electrons with eta < DET->MAX_CALO_FWD
-                DET->SmearElectron(genMomentum);
+                nonS = genMomentum.E();
+                DET->SmearElectron(genMomentum);
+                if(eta < DET->MAX_TRACKER){
+                elementElec=(RESOLELEC*) branchelec->NewEntry();
+                elementElec->NonSmeareE = nonS;
+                elementElec->SmeareE = genMomentum.E();}
                 break; // case pE
               case pGAMMA: // all photons with eta < DET->MAX_CALO_FWD
@@ -251 +315 @@
                 break; // case pGAMMA
               case pMU: // all muons with eta < DET->MAX_MU
+                elementMuon = (RESOLMUON*) branchmuon->NewEntry();
+                elementMuon->OneoverNonSmearePT = (1/genMomentum.Pt());
                 DET->SmearMu(genMomentum);
+                elementMuon->OneoverSmearePT = (1/genMomentum.Pt());
                 break; // case pMU
               case pLAMBDA: // all lambdas with eta < DET->MAX_CALO_FWD
@@ -264 +331 @@
               if(pid != pMU)
                 {
-//                PTmisS = PTmisS + genMomentum;
                   towers.push_back(PhysicsTower(LorentzVector(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E())));
-                  if(genMomentum.Et()>0.5)input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
+                  input_particlesS.push_back(fastjet::PseudoJet(genMomentum.Px(),genMomentum.Py(),genMomentum.Pz(), genMomentum.E()));
                 }
               
@@ -286 +352 @@
             } // switch
         } // while
-     
-     TLorentzVector TowerTLV(0.,0.,0.,0.);
+    
+     TLorentzVector Att(0.,0.,0.,0.);
+     float ScalarEt=0;
      for(unsigned int i=0; i < towers.size(); i++) 
        {
-         TowerTLV.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
-         if(TowerTLV.Et()>0.5)PTmisS = PTmisS + TowerTLV;
-      }
-
- 
+         Att.SetPxPyPzE(towers[i].fourVector.px,towers[i].fourVector.py,towers[i].fourVector.pz,towers[i].fourVector.E);
+         ScalarEt = ScalarEt + Att.Et();
+         PTmisS = PTmisS + Att;
+       }
+//cout<<"non smeare "<<PTmis.Pt()<<" "<<PTmisS.Pt()<<endl;
+//cout<<"smeare "<<PTmis.Px()<<" "<<PTmisS.Px()<<endl;
+//cout<<"**********"<<endl;
+
       elementEtmis= (ETMIS*) branchetmis->NewEntry();
-      elementEtmis->Et = (PTmis).Et();
-      elementEtmis->EtSmeare = (PTmisS).Et()-(PTmis).Et();
+      elementEtmis->Et = (PTmis).Pt();
+      elementEtmis->Ex = (-PTmis).Px();
+      elementEtmis->SEt = ScalarEt;
+
+      elementEtmis->EtSmeare = (PTmisS).Pt()-(PTmis).Pt();
+      elementEtmis->ExSmeare = (-PTmisS).Px()-(PTmis).Px();
       
       //*****************************
@@ -316 +390 @@
         }
       
+      TSimpleArray<TRootGenParticle> TausHadr = TauHadr(branchGen);
+
       TLorentzVector JETSm(0,0,0,0);
       for (unsigned int i = 0; i < sorted_jets.size(); i++) {
@@ -326 +402 @@
             elementJet->NonSmearePT = JET.Et();
             elementJet->SmearePT = JETSm.Et()/JET.Et();
-            
           }
+      }
+      for (unsigned int i = 0; i < sorted_jetsS.size(); i++) {
+        TLorentzVector JETT(0,0,0,0);
+        JETT.SetPxPyPzE(sorted_jetsS[i].px(),sorted_jetsS[i].py(),sorted_jetsS[i].pz(),sorted_jetsS[i].E());
+        if(fabs(JETT.Eta()) < (DET->MAX_TRACKER - DET->TAU_CONE_TRACKS)) 
+           {
+            for(Int_t i=0; i<TausHadr.GetEntries();i++)
+              {
+                if(DeltaR(TausHadr[i]->Phi,TausHadr[i]->Eta,JETT.Phi(),JETT.Eta())<0.1)
+                  {
+                    elementTaujet= (TAUHAD*) branchtaujet->NewEntry();
+                    elementTaujet->EnergieCen = (DET->EnergySmallCone(towers,JETT.Eta(),JETT.Phi())/JETT.E());
+                    elementTaujet->NumTrack = DET->NumTracks(TrackCentral,DET->PT_TRACK_TAU,JETT.Eta(),JETT.Phi());
+                 }
+              }
+           }
+         
         
       } // for itJet : loop on all jets