Added file check.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderESD.cxx
index 68f0568aabc95dc9850abe906f642a0d5d26b5c0..d58d34a657f04651fd8c9c121a90ee98c60d5816 100644 (file)
@@ -11,6 +11,7 @@
 //                                                                  //
 //////////////////////////////////////////////////////////////////////
 
+#include <Riostream.h>
 #include <TMath.h>
 #include <TPDGCode.h>
 #include <TString.h>
 #include <TH1.h>
 #include <AliESD.h>
 #include <AliESDtrack.h>
+#include <AliKalmanTrack.h>
 #include <AliJetEventParticles.h>
 #include "AliJetParticlesReaderESD.h"
 
 ClassImp(AliJetParticlesReaderESD)
 
-AliJetParticlesReaderESD::AliJetParticlesReaderESD(const Char_t* esdfilename) :
+AliJetParticlesReaderESD::AliJetParticlesReaderESD(Bool_t constrained,
+                                                  const Char_t* esdfilename) :
   AliJetParticlesReader(),
+  fConstrained(constrained),
   fESDFileName(esdfilename),
+  fESD(0),
   fFile(0),
+  fTree(0),
   fKeyIterator(0),
   fPassFlag(AliESDtrack::kTPCrefit)
 {
@@ -39,11 +45,15 @@ AliJetParticlesReaderESD::AliJetParticlesReaderESD(const Char_t* esdfilename) :
 /********************************************************************/
   
 AliJetParticlesReaderESD::AliJetParticlesReaderESD(
+                                      Bool_t constrained,
                                       TObjArray* dirs,
                                       const Char_t* esdfilename) :
   AliJetParticlesReader(dirs),
+  fConstrained(constrained),
   fESDFileName(esdfilename),
+  fESD(0),
   fFile(0),
+  fTree(0),
   fKeyIterator(0),
   fPassFlag(AliESDtrack::kTPCrefit)
 {
@@ -55,69 +65,10 @@ AliJetParticlesReaderESD::AliJetParticlesReaderESD(
 AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
 {
   //desctructor
-  if(fFile) delete fFile;
+  if(fESD) delete fESD;
+  if(fTree) delete fTree;
   if(fKeyIterator) delete fKeyIterator;
-}
-
-
-Int_t AliJetParticlesReaderESD::ReadNext()
-{
-  //reads next event from fFile
-
-  do   // is OK even if 0 dirs specified, 
-    {  // in that case we try to read from "./"
-      if (fFile == 0)
-       {
-         fFile = OpenFile(fCurrentDir);
-         if (fFile == 0)
-           {
-             Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
-             fCurrentDir++;
-             continue;
-           }
-     
-         fKeyIterator = new TIter(fFile->GetListOfKeys());  
-         fCurrentEvent = 0;
-         //fFile->Dump();
-         //fFile->GetListOfKeys()->Print();
-       } 
-
-      TKey* key = (TKey*)fKeyIterator->Next();
-      if (key == 0)
-       {
-         fCurrentDir++;
-         delete fKeyIterator;
-         fKeyIterator = 0;
-         delete fFile; //we have to assume there are no more ESD objects in the fFile
-         fFile = 0;
-         continue;
-       }
-
-      TString esdname = "ESD";
-      esdname+=fCurrentEvent;
-      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
-      if (esd == 0)
-      {
-       //Info("ReadNext","This key is not an AliESD object %s",key->GetName());
-       Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
-       
-       fCurrentDir++;
-       delete fKeyIterator;
-       fKeyIterator = 0;
-       delete fFile;//we have to assume there is no more ESD objects in the fFile
-       fFile = 0;
-       continue;
-      }
-     
-      ReadESD(esd);
-      
-      fCurrentEvent++;
-      fNEventsRead++;
-      delete esd;
-      return kTRUE;//success -> read one event
-    }  while(fCurrentDir < GetNumberOfDirs());
-       //end of loop over directories specified in fDirs Obj Array  
-  return kFALSE; //no more directories to read
+  if(fFile) delete fFile;
 }
 
 /**********************************************************/
@@ -139,17 +90,29 @@ Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
   //   return kFALSE;
   // }
    
-  //Float_t mf = esd->GetMagneticField(); 
-  //if ( (mf == 0.0) && (fNTrackPoints > 0) )
-  //{
-  //   Error("ReadESD","Magnetic Field is 0 and Track Points demanded. Skipping to next event.");
-  //   return kFALSE;
-  //}
+  Float_t mf = esd->GetMagneticField(); 
+  if (mf <= 0.0)  
+  {
+     Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
+     return kFALSE;
+  }
+  AliKalmanTrack::SetMagneticField(mf/10.);
 
   Info("ReadESD","Reading Event %d",fCurrentEvent);
   if((!fOwner) || (fEventParticles==0)) 
     fEventParticles = new AliJetEventParticles();
 
+  const Int_t kntr = esd->GetNumberOfTracks();
+  Info("ReadESD","Found %d tracks.",kntr);
+  fEventParticles->Reset(kntr);
+
+  TString headdesc="";
+  headdesc+="Run ";
+  headdesc+=esd->GetRunNumber();
+  headdesc+=": Ev ";
+  headdesc+=esd->GetEventNumber();
+  fEventParticles->SetHeader(headdesc);
+
   Double_t vertexpos[3];//vertex position
   const AliESDVertex* kvertex = esd->GetVertex();
   if (kvertex == 0)
@@ -165,29 +128,33 @@ Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
    }
   fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
 
-  const Int_t kntr = esd->GetNumberOfTracks();
-  Info("ReadESD","Found %d tracks.",kntr);
+  //loop over tracks
   for (Int_t i = 0;i<kntr; i++)
    {
+
      const AliESDtrack *kesdtrack = esd->GetTrack(i);
      if (kesdtrack == 0)
       {
-        Error("ReadESD","Can not get track %d", i);
+        Error("ReadESD","Can't get track %d", i);
         continue;
       }
 
-     if ((kesdtrack->GetStatus() & fPassFlag) == kFALSE)
+     if ((kesdtrack->GetStatus() & fPassFlag)!=fPassFlag)
       {
-       Info("ReadNext","Particle skipped.");
+       Info("ReadNext","Particle skipped: %u.",kesdtrack->GetStatus());
         continue;
       }
 
      Double_t mom[3];  //momentum
-     kesdtrack->GetPxPyPz(mom);
-     //kesdtrack->GetConstrainedPxPyPz(mom);
-     //Double_t pos[3];//position
-     //kesdtrack->GetXYZ(pos);
-     //kesdtrack->GetConstrainedXYZ(pos);
+     Double_t xyz[3];  //position
+     if (fConstrained) {
+       if (kesdtrack->GetConstrainedChi2() > 25) continue;
+       kesdtrack->GetConstrainedPxPyPz(mom);
+       kesdtrack->GetConstrainedXYZ(xyz);
+     } else {
+       kesdtrack->GetPxPyPz(mom);
+       kesdtrack->GetXYZ(xyz);
+     }
      const Float_t kmass=kesdtrack->GetMass();
      const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
      const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
@@ -195,8 +162,16 @@ Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
      const Float_t kp=TMath::Sqrt(kp2);
      const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30)); 
      const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
+     //Double_t dx = xyz[0]-vertexpos[0];
+     //Double_t dy = xyz[1]-vertexpos[1];
+     //Float_t dca = TMath::Sqrt(dx*dx + dy*dy);
+     //Float_t dz = xyz[2]-vertexpos[2];
+     UInt_t index[6];
+     const Int_t kncl=kesdtrack->GetITSclusters(index)
+                      +kesdtrack->GetTPCclusters(NULL)
+                      +kesdtrack->GetTRDclusters(NULL);
      if(IsAcceptedParticle(kpt,kphi,keta))
-       fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kpt,kphi,keta);
+       fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kncl,kpt,kphi,keta);
 
    } // loop over tracks
 
@@ -219,6 +194,86 @@ void AliJetParticlesReaderESD::Rewind()
 
 /**********************************************************/
 
+Int_t AliJetParticlesReaderESD::ReadNext()
+{
+  //reads next event from fFile
+
+  do   // is OK even if 0 dirs specified, 
+    {  // in that case we try to read from "./"
+      if (fFile == 0)
+       {
+         fFile = OpenFile(fCurrentDir);
+         if (fFile == 0)
+           {
+             Error("ReadNext","Can't get fFile for dir no. %d",fCurrentDir);
+             fCurrentDir++;
+             continue;
+           }
+     
+         fCurrentEvent = 0;
+         //fFile->GetListOfKeys()->Print();
+         
+         if(fTree) delete fTree;
+         fTree = dynamic_cast<TTree*>(fFile->Get("esdTree"));
+         if(fTree)
+           fTree->SetBranchAddress("ESD",&fESD);
+         else
+           fKeyIterator = new TIter(fFile->GetListOfKeys());  
+       } 
+
+      if(fTree)
+       {
+         if(AliKalmanTrack::GetConvConst()<=0.)
+           AliKalmanTrack::SetMagneticField(0.5);
+         if(fCurrentEvent>=fTree->GetEntries())
+           {
+             fCurrentDir++;
+             delete fTree;
+             fTree = 0;
+             delete fFile;
+             fFile = 0;
+             continue;
+           }
+         fTree->GetEvent(fCurrentEvent);
+       } 
+      else 
+       { // "old" way via ESD objects stored in root file
+         TKey* key = (TKey*)fKeyIterator->Next();
+         if (key == 0)
+           {
+             fCurrentDir++;
+             delete fKeyIterator;
+             fKeyIterator = 0;
+             delete fFile; //we have to assume there are no more ESD objects in the fFile
+             fFile = 0;
+             continue;
+           }
+         TString esdname = "ESD";
+         esdname+=fCurrentEvent;
+         if(fESD) delete fESD;
+         fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
+         if (fESD == 0)
+           {
+             Info("ReadNext","Can't find AliESD object named %s",esdname.Data());
+             fCurrentDir++;
+             delete fKeyIterator;
+             fKeyIterator = 0;
+             delete fFile;//we have to assume there is no more ESD objects in the fFile
+             fFile = 0;
+             continue;
+           }
+       }
+      ReadESD(fESD);
+      fCurrentEvent++;
+      fNEventsRead++;
+      return kTRUE;//success -> read one event
+    }  while(fCurrentDir < GetNumberOfDirs());
+      //end of loop over directories specified in fDirs Obj Array  
+  return kFALSE; //no more directories to read
+}
+
+/**********************************************************/
+
 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
 {
   //opens fFile with kine tree
@@ -226,7 +281,7 @@ TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
  const TString& dirname = GetDirName(n);
  if (dirname == "")
   {
-   Error("OpenFiles","Can not get directory name");
+   Error("OpenFiles","Can't get directory name");
    return 0;
   }
  TString filename = dirname +"/"+ fESDFileName;