]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ANALYSIS/AliReaderESD.cxx
Additional safety check introduced
[u/mrichter/AliRoot.git] / ANALYSIS / AliReaderESD.cxx
index 4d2e75bea8a174e47c2972698fb0dad942b5099a..17e09a737d065511382c807bfc2a584a3a7cd28f 100644 (file)
 #include <AliRunLoader.h>
 #include <AliStack.h>
 #include <AliESDtrack.h>
+#include <AliKalmanTrack.h>
 #include <AliESD.h>
 
 #include "AliAnalysis.h"
 #include "AliAODRun.h"
 #include "AliAOD.h"
-#include "AliAODStdParticle.h"
+#include "AliAODParticle.h"
 #include "AliAODParticleCut.h"
 #include "AliTrackPoints.h"
 #include "AliClusterMap.h"
@@ -47,6 +48,8 @@ AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename)
  fNTrackPoints(0),
  fdR(0.0),
  fClusterMap(kFALSE),
+ fITSTrackPoints(kFALSE),
+ fMustTPC(kFALSE),
  fNTPCClustMin(0),
  fNTPCClustMax(150),
  fTPCChi2PerClustMin(0.0),
@@ -94,6 +97,8 @@ AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char
  fNTrackPoints(0),
  fdR(0.0),
  fClusterMap(kFALSE),
+ fITSTrackPoints(kFALSE),
+ fMustTPC(kFALSE),
  fNTPCClustMin(0),
  fNTPCClustMax(150),
  fTPCChi2PerClustMin(0.0),
@@ -140,7 +145,7 @@ Int_t AliReaderESD::ReadNext()
 //reads next event from fFile
   //fRunLoader is for reading Kine
   
-  if (AliAODParticle::GetDebug())
+  if (AliVAODParticle::GetDebug())
     Info("ReadNext","Entered");
     
   if (fEventSim == 0x0)  fEventSim = new AliAOD();
@@ -168,7 +173,7 @@ Int_t AliReaderESD::ReadNext()
      TKey* key = (TKey*)fKeyIterator->Next();
      if (key == 0x0)
       {
-        if (AliAODParticle::GetDebug() > 2 )
+        if (AliVAODParticle::GetDebug() > 2 )
           {
             Info("ReadNext","No more keys.");
           }
@@ -187,7 +192,7 @@ Int_t AliReaderESD::ReadNext()
 //     TObject* esdobj = key->ReadObj();
 //     if (esdobj == 0x0)
 //      {
-//        if (AliAODParticle::GetDebug() > 2 )
+//        if (AliVAODParticle::GetDebug() > 2 )
 //          {
 //            Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
 //            key->Dump();
@@ -202,11 +207,11 @@ Int_t AliReaderESD::ReadNext()
      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
      if (esd == 0x0)
       {
-//        if (AliAODParticle::GetDebug() > 2 )
+//        if (AliVAODParticle::GetDebug() > 2 )
 //          {
 //            Info("ReadNext","This key is not an AliESD object %s",key->GetName());
 //          }
-        if (AliAODParticle::GetDebug() > 2 )
+        if (AliVAODParticle::GetDebug() > 2 )
           {
             Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
           }
@@ -248,7 +253,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
      Error("ReadESD","ESD is NULL");
      return 1;
    }
-
+  
   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
   if (pdgdb == 0x0)
    {
@@ -258,12 +263,18 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
    
   Float_t mf = esd->GetMagneticField(); 
 
-  if ( (mf == 0.0) && (fNTrackPoints > 0) )
+  if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
    {
       Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
       return 1;
    }
 
+  if (fITSTrackPoints)
+   {
+     Info("ReadESD","Magnetic Field is %f",mf/10.);
+     AliKalmanTrack::SetMagneticField(mf/10.);
+   }
   AliStack* stack = 0x0;
   if (fReadSim && fRunLoader)
    {
@@ -284,7 +295,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
      vertex->GetXYZ(vertexpos);
    }
    
-  if (AliAODParticle::GetDebug() > 0)
+  if (AliVAODParticle::GetDebug() > 0)
    {
      Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]);
    }
@@ -305,14 +316,24 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
      //if (esdtrack->HasVertexParameters() == kFALSE) 
      if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
       {
-        if (AliAODParticle::GetDebug() > 2) 
+        if (AliVAODParticle::GetDebug() > 2) 
           Info("ReadNext","Particle skipped: Data at vertex not available.");
         continue;
       }
 
+     if (fMustTPC)
+      {
+       if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
+        {
+          if (AliVAODParticle::GetDebug() > 2) 
+            Info("ReadNext","Particle skipped: Was not reconstructed in TPC.");
+          continue;
+        }
+      }     
+
      if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE) 
       {
-        if (AliAODParticle::GetDebug() > 2) 
+        if (AliVAODParticle::GetDebug() > 2) 
           Info("ReadNext","Particle skipped: PID BIT is not set.");
         continue;
       }
@@ -323,7 +344,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
      esdtrack->GetConstrainedExternalParameters(extx,extp);
      if (extp[4] == 0.0)
       {
-        if (AliAODParticle::GetDebug() > 2) 
+        if (AliVAODParticle::GetDebug() > 2) 
           Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
         continue;
       } 
@@ -337,7 +358,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
      Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
 
      //Particle from kinematics
-     AliAODStdParticle* particle = 0;
+     AliAODParticle* particle = 0;
      Bool_t keeppart = kFALSE;
      if ( fReadSim && stack  )
       {
@@ -348,14 +369,18 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
            Error("ReadNext","Can not find track with such label.");
            continue;
          }
-        if(Pass(p->GetPdgCode())) 
+         
+        if (fCheckParticlePID)
          {
-           if ( AliAODParticle::GetDebug() > 5 )
-             Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode());
-           continue; //check if we are intersted with particles of this type 
-         }
+           if(Rejected(p->GetPdgCode())) 
+            {
+              if ( AliVAODParticle::GetDebug() > 5 )
+                Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode());
+              continue; //check if we are intersted with particles of this type 
+            }
+         }  
 //           if(p->GetPdgCode()<0) charge = -1;
-        particle = new AliAODStdParticle(*p,i);
+        particle = new AliAODParticle(*p,i);
 
       }
       
@@ -366,14 +391,14 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
      for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
      if (rc==0.0) 
       {
-        if (AliAODParticle::GetDebug() > 2) 
+        if (AliVAODParticle::GetDebug() > 2) 
           Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
         continue;
       }
 
      for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
 
-     if (AliAODParticle::GetDebug() > 4)
+     if (AliVAODParticle::GetDebug() > 4)
       { 
         Info("ReadNext","###########################################################################");
         Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
@@ -392,12 +417,20 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
            msg+=")";
          }
         Info("ReadNext","%s",msg.Data());
-      }//if (AliAODParticle::GetDebug()>4)
+      }//if (AliVAODParticle::GetDebug()>4)
 
       AliTrackPoints* tpts = 0x0;
       if (fNTrackPoints > 0) 
        {
          tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
+//         tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
+       }
+
+      AliTrackPoints* itstpts = 0x0;
+      if (fITSTrackPoints) 
+       {
+         itstpts = new AliTrackPoints(AliTrackPoints::kITS,esdtrack);
+//         itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
        }
 
       AliClusterMap* cmap = 0x0; 
@@ -437,14 +470,14 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
         Float_t pp = w[s];
         if (pp == 0.0) 
          {
-           if ( AliAODParticle::GetDebug() > 5 )
+           if ( AliVAODParticle::GetDebug() > 5 )
              Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
            continue;
          }
 
-        if(Pass(pdgcode)) 
+        if(Rejected(pdgcode)) 
          {
-           if ( AliAODParticle::GetDebug() > 5 )
+           if ( AliVAODParticle::GetDebug() > 5 )
              Info("ReadNext","PID (%d) did not pass the cut.",pdgcode);
            continue; //check if we are intersted with particles of this type 
          }
@@ -452,7 +485,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
         Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
         Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
 
-        AliAODStdParticle* track = new AliAODStdParticle(pdgcode, w[s],i, 
+        AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i, 
                                                    mom[0], mom[1], mom[2], tEtot,
                                                    pos[0], pos[1], pos[2], 0.);
         //copy probabilitis of other species (if not zero)
@@ -463,10 +496,10 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
            track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
          }
 
-        if(Pass(track))//check if meets all criteria of any of our cuts
+        if(Rejected(track))//check if meets all criteria of any of our cuts
                        //if it does not delete it and take next good track
          { 
-           if ( AliAODParticle::GetDebug() > 4 )
+           if ( AliVAODParticle::GetDebug() > 4 )
              Info("ReadNext","Track did not pass the cut");
            delete track;
            continue;
@@ -475,7 +508,12 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
         //Single Particle cuts on cluster map and track points rather do not have sense
         if (tpts)
          {
-           track->SetTrackPoints(tpts); 
+           track->SetTPCTrackPoints(tpts);
+         }
+         
+        if (itstpts)
+         {
+           track->SetITSTrackPoints(itstpts); 
          }
 
         if (cmap) 
@@ -487,7 +525,7 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
         if (particle) fEventSim->AddParticle(particle);
         keeppart = kTRUE;
 
-        if (AliAODParticle::GetDebug() > 4 )
+        if (AliVAODParticle::GetDebug() > 4 )
          {
            Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
            track->Print();
@@ -500,8 +538,26 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
       {
         delete particle;//particle was not stored in event
         delete tpts;
+        delete itstpts;
         delete cmap;
       }
+     else
+      {
+        if ( fReadSim && stack  )
+         {
+           if (particle->P() < 0.00001)
+            {
+              Info("ReadNext","###################################");
+              Info("ReadNext","###################################");
+              Info("ReadNext","Track Label %d",esdtrack->GetLabel());
+              TParticle *p = stack->Particle(esdtrack->GetLabel());
+              Info("ReadNext","");
+              p->Print();
+              Info("ReadNext","");
+              particle->Print();
+            }
+         }   
+      } 
 
    }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
 
@@ -509,6 +565,15 @@ Int_t AliReaderESD::ReadESD(AliESD* esd)
          fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(),
          fNEventsRead,fCurrentEvent,fCurrentDir);
   fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
+  
+  /******************************************************/
+  /******     Setting glevet properties     *************/
+  /******************************************************/
+    
+  if (fEventRec->GetNumberOfParticles() > 0)
+   {
+     fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
+   }
   return 0;
 }