Bugfixes again for HLT.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderESD.cxx
index dbf28dd8bfde0762a8c8bc0102fbe5ffa433388b..735e174364a73514370835a946968cfa6a6b911d 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),
@@ -41,9 +45,11 @@ 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),
@@ -59,12 +65,138 @@ AliJetParticlesReaderESD::AliJetParticlesReaderESD(
 AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
 {
   //desctructor
-  if(fFile) delete fFile;
   if(fESD) delete fESD;
   if(fTree) delete fTree;
   if(fKeyIterator) delete fKeyIterator;
+  if(fFile) delete fFile;
+}
+
+/**********************************************************/
+
+Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
+{
+  //Reads one ESD
+
+  if (esd == 0)
+   {
+     Error("ReadESD","ESD is NULL");
+     return kFALSE;
+   }
+
+  /*
+  TDatabasePDG* pdgdb = TDatabasePDG::Instance();
+  if (pdgdb == 0)
+  {
+     Error("ReadESD","Can not get PDG Database Instance.");
+     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",fCurrentDir*1000+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)
+   {
+     Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
+     vertexpos[0] = 0.0;
+     vertexpos[1] = 0.0;
+     vertexpos[2] = 0.0;
+   }
+  else
+   {
+     kvertex->GetXYZ(vertexpos);
+   }
+  fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
+
+  //loop over tracks
+  for (Int_t i = 0;i<kntr; i++)
+   {
+
+     const AliESDtrack *kesdtrack = esd->GetTrack(i);
+     if (kesdtrack == 0)
+      {
+        Error("ReadESD","Can't get track %d", i);
+        continue;
+      }
+     
+     ULong_t cmptest=(kesdtrack->GetStatus() & fPassFlag);
+     if (cmptest!=fPassFlag)
+      {
+       Info("ReadNext","Particle %d skipped: %u.",i,kesdtrack->GetStatus());
+       cout << i << " "; PrintESDtrack(kesdtrack); cout << endl;
+        continue;
+      }
+
+     Double_t mom[3];  //momentum
+     Double_t xyz[3];  //position
+     if (fConstrained) {
+       if (kesdtrack->GetConstrainedChi2() > 25) continue;
+       kesdtrack->GetConstrainedPxPyPz(mom);
+       kesdtrack->GetConstrainedXYZ(xyz);
+     } else {
+       if(!kesdtrack->GetPxPyPzAt(0,mom)) continue;
+       kesdtrack->GetXYZAt(0, 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);
+     const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+     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(),kncl,kpt,kphi,keta);
+
+   } // loop over tracks
+
+  return kTRUE;
+}
+
+/**********************************************************/
+
+void AliJetParticlesReaderESD::Rewind()
+{
+  //rewinds reading 
+
+  if(fFile) delete fFile;
+  if(fKeyIterator) delete fKeyIterator;
+  fKeyIterator = 0;
+  fFile = 0;
+  fCurrentDir = 0;
+  fNEventsRead = 0;
 }
 
+/**********************************************************/
 
 Int_t AliJetParticlesReaderESD::ReadNext()
 {
@@ -77,7 +209,7 @@ Int_t AliJetParticlesReaderESD::ReadNext()
          fFile = OpenFile(fCurrentDir);
          if (fFile == 0)
            {
-             Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
+             Error("ReadNext","Can't get fFile for dir no. %d",fCurrentDir);
              fCurrentDir++;
              continue;
            }
@@ -95,6 +227,8 @@ Int_t AliJetParticlesReaderESD::ReadNext()
 
       if(fTree)
        {
+         if(AliKalmanTrack::GetConvConst()<=0.)
+           AliKalmanTrack::SetMagneticField(0.5);
          if(fCurrentEvent>=fTree->GetEntries())
            {
              fCurrentDir++;
@@ -124,7 +258,7 @@ Int_t AliJetParticlesReaderESD::ReadNext()
          fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
          if (fESD == 0)
            {
-             Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
+             Info("ReadNext","Can't find AliESD object named %s",esdname.Data());
              fCurrentDir++;
              delete fKeyIterator;
              fKeyIterator = 0;
@@ -144,103 +278,6 @@ Int_t AliJetParticlesReaderESD::ReadNext()
 
 /**********************************************************/
 
-Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
-{
-  //Reads one ESD
-
-  if (esd == 0)
-   {
-     Error("ReadESD","ESD is NULL");
-     return kFALSE;
-   }
-
-  //TDatabasePDG* pdgdb = TDatabasePDG::Instance();
-  //if (pdgdb == 0)
-  //{
-  //   Error("ReadESD","Can not get PDG Database Instance.");
-  //   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;
-  //}
-
-  Info("ReadESD","Reading Event %d",fCurrentEvent);
-  if((!fOwner) || (fEventParticles==0)) 
-    fEventParticles = new AliJetEventParticles();
-
-  Double_t vertexpos[3];//vertex position
-  const AliESDVertex* kvertex = esd->GetVertex();
-  if (kvertex == 0)
-   {
-     Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
-     vertexpos[0] = 0.0;
-     vertexpos[1] = 0.0;
-     vertexpos[2] = 0.0;
-   }
-  else
-   {
-     kvertex->GetXYZ(vertexpos);
-   }
-  fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
-
-  const Int_t kntr = esd->GetNumberOfTracks();
-  Info("ReadESD","Found %d tracks.",kntr);
-  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);
-        continue;
-      }
-
-     if ((kesdtrack->GetStatus() & fPassFlag) == kFALSE)
-      {
-       Info("ReadNext","Particle skipped: %ud.",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);
-     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);
-     const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
-     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]);
-     if(IsAcceptedParticle(kpt,kphi,keta))
-       fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kpt,kphi,keta);
-
-   } // loop over tracks
-
-  return kTRUE;
-}
-
-/**********************************************************/
-
-void AliJetParticlesReaderESD::Rewind()
-{
-  //rewinds reading 
-
-  if(fFile) delete fFile;
-  if(fKeyIterator) delete fKeyIterator;
-  fKeyIterator = 0;
-  fFile = 0;
-  fCurrentDir = 0;
-  fNEventsRead = 0;
-}
-
-/**********************************************************/
-
 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
 {
   //opens fFile with kine tree
@@ -248,7 +285,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;
@@ -266,3 +303,45 @@ TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
    
  return ret;
 }
+
+/**********************************************************/
+
+void AliJetParticlesReaderESD::PrintESDtrack(const AliESDtrack *kesdtrack) const
+{
+  ULong_t status=kesdtrack->GetStatus();
+  cout << hex << status << dec << ": ";
+  if((status & AliESDtrack::kITSin) == AliESDtrack::kITSin) cout << "ITSin ";
+  if((status & AliESDtrack::kITSout) == AliESDtrack::kITSout) cout << "ITSout ";
+  if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) cout << "ITSrefit ";
+  if((status & AliESDtrack::kITSpid) == AliESDtrack::kITSpid) cout << "ITSpid ";
+
+  if((status & AliESDtrack::kTPCin)  == AliESDtrack::kTPCin) cout << "TPCin ";
+  if((status & AliESDtrack::kTPCout) == AliESDtrack::kTPCout) cout << "TPCout ";
+  if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) cout << "TPCrefit ";
+  if((status & AliESDtrack::kTPCpid) == AliESDtrack::kTPCpid) cout << "TPCpid ";
+
+  if((status & AliESDtrack::kTRDin) == AliESDtrack::kTRDin) cout << "TRDin ";
+  if((status & AliESDtrack::kTRDout) == AliESDtrack::kTRDout) cout << "TRDout ";
+  if((status & AliESDtrack::kTRDrefit) == AliESDtrack::kTRDrefit) cout << "TRDrefit ";
+  if((status & AliESDtrack::kTRDpid) == AliESDtrack::kTRDpid) cout << "TRDpid ";
+  if((status & AliESDtrack::kTRDbackup) == AliESDtrack::kTRDbackup) cout << "TRDbackup ";
+  if((status & AliESDtrack::kTRDStop) == AliESDtrack::kTRDStop) cout << "TRDstop ";
+
+  if((status & AliESDtrack::kTOFin) == AliESDtrack::kTOFin) cout << "TOFin ";
+  if((status & AliESDtrack::kTOFout) == AliESDtrack::kTOFout) cout << "TOFout ";
+  if((status & AliESDtrack::kTOFrefit) == AliESDtrack::kTOFrefit) cout << "TOFrefit ";
+  if((status & AliESDtrack::kTOFpid) == AliESDtrack::kTOFpid) cout << "TOFpid ";
+
+  if((status & AliESDtrack::kPHOSpid) == AliESDtrack::kPHOSpid) cout << "PHOSpid ";
+  if((status & AliESDtrack::kRICHpid) == AliESDtrack::kRICHpid) cout << "RICHpid ";
+  if((status & AliESDtrack::kEMCALpid) == AliESDtrack::kEMCALpid) cout << "EMCALpid ";
+  if((status & AliESDtrack::kESDpid) == AliESDtrack::kESDpid) cout << "ESDpid ";
+  if((status & AliESDtrack::kTIME) == AliESDtrack::kTIME) cout << "TIME ";
+  cout << endl; 
+  /*
+  cout << kesdtrack->GetConstrainedChi2() 
+       << " " << kesdtrack->GetITSchi2()
+       << " " << kesdtrack->GetTPCchi2() 
+       << " " << kesdtrack->GetTRDchi2()<< endl;
+  */
+}