JETAN/AliJetParticlesReaderESD.cxx

   1 // $Id$
   2
   3 //____________________________________________________________________
   4 //////////////////////////////////////////////////////////////////////
   5 //                                                                  //
   6 // class AliHBTReaderESD                                            //
   7 //                                                                  //
   8 // Reader for ALICE Event Summary Data (ESD).                       //
   9 //                                                                  //
  10 // Piotr.Skowronski@cern.ch                                         //
  11 //                                                                  //
  12 //////////////////////////////////////////////////////////////////////
  13
  14 #include <Riostream.h>
  15 #include <TMath.h>
  16 #include <TPDGCode.h>
  17 #include <TString.h>
  18 #include <TObjString.h>
  19 #include <TTree.h>
  20 #include <TFile.h>
  21 #include <TKey.h>
  22 #include <TH1.h>
  23 #include <AliESD.h>
  24 #include <AliESDtrack.h>
  25 #include <AliKalmanTrack.h>
  26 #include <AliJetEventParticles.h>
  27 #include "AliJetParticlesReaderESD.h"
  28
  29 ClassImp(AliJetParticlesReaderESD)
  30
  31 AliJetParticlesReaderESD::AliJetParticlesReaderESD(Bool_t constrained,
  32                                                    const Char_t* esdfilename) :
  33   AliJetParticlesReader(),
  34   fConstrained(constrained),
  35   fESDFileName(esdfilename),
  36   fESD(0),
  37   fFile(0),
  38   fTree(0),
  39   fKeyIterator(0),
  40   fPassFlag(AliESDtrack::kTPCrefit)
  41 {
  42   //constructor
  43 }
  44
  45 /********************************************************************/
  46
  47 AliJetParticlesReaderESD::AliJetParticlesReaderESD(
  48                                       Bool_t constrained,
  49                                       TObjArray* dirs,
  50                                       const Char_t* esdfilename) :
  51   AliJetParticlesReader(dirs),
  52   fConstrained(constrained),
  53   fESDFileName(esdfilename),
  54   fESD(0),
  55   fFile(0),
  56   fTree(0),
  57   fKeyIterator(0),
  58   fPassFlag(AliESDtrack::kTPCrefit)
  59 {
  60   //constructor
  61 }
  62
  63 /********************************************************************/
  64
  65 AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
  66 {
  67   //desctructor
  68   if(fESD) delete fESD;
  69   if(fTree) delete fTree;
  70   if(fKeyIterator) delete fKeyIterator;
  71   if(fFile) delete fFile;
  72 }
  73
  74 /**********************************************************/
  75
  76 Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
  77 {
  78   //Reads one ESD
  79
  80   if (esd == 0)
  81    {
  82      Error("ReadESD","ESD is NULL");
  83      return kFALSE;
  84    }
  85
  86   /*
  87   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
  88   if (pdgdb == 0)
  89   {
  90      Error("ReadESD","Can not get PDG Database Instance.");
  91      return kFALSE;
  92   }
  93   */
  94
  95   Float_t mf = esd->GetMagneticField();
  96   if (mf <= 0.0)
  97   {
  98      Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
  99      return kFALSE;
 100   }
 101   AliKalmanTrack::SetMagneticField(mf/10.);
 102
 103   Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
 104   if((!fOwner) || (fEventParticles==0))
 105     fEventParticles = new AliJetEventParticles();
 106
 107   const Int_t kntr = esd->GetNumberOfTracks();
 108   Info("ReadESD","Found %d tracks.",kntr);
 109   fEventParticles->Reset(kntr);
 110
 111   TString headdesc="";
 112   headdesc+="Run ";
 113   headdesc+=esd->GetRunNumber();
 114   headdesc+=": Ev ";
 115   headdesc+=esd->GetEventNumber();
 116   fEventParticles->SetHeader(headdesc);
 117
 118   Double_t vertexpos[3];//vertex position
 119   const AliESDVertex* kvertex = esd->GetVertex();
 120   if (kvertex == 0)
 121    {
 122      Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
 123      vertexpos[0] = 0.0;
 124      vertexpos[1] = 0.0;
 125      vertexpos[2] = 0.0;
 126    }
 127   else
 128    {
 129      kvertex->GetXYZ(vertexpos);
 130    }
 131   fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
 132
 133   //loop over tracks
 134   for (Int_t i = 0;i<kntr; i++)
 135    {
 136
 137      const AliESDtrack *kesdtrack = esd->GetTrack(i);
 138      if (kesdtrack == 0)
 139       {
 140         Error("ReadESD","Can't get track %d", i);
 141         continue;
 142       }
 143
 144      ULong_t cmptest=(kesdtrack->GetStatus() & fPassFlag);
 145      if (cmptest!=fPassFlag)
 146       {
 147         Info("ReadNext","Particle %d skipped: %u.",i,kesdtrack->GetStatus());
 148         cout << i << " "; PrintESDtrack(kesdtrack); cout << endl;
 149         continue;
 150       }
 151
 152      Double_t mom[3];  //momentum
 153      Double_t xyz[3];  //position
 154      if (fConstrained) {
 155        if (kesdtrack->GetConstrainedChi2() > 25) continue;
 156        kesdtrack->GetConstrainedPxPyPz(mom);
 157        kesdtrack->GetConstrainedXYZ(xyz);
 158      } else {
 159        if(!kesdtrack->GetPxPyPzAt(0,mom)) continue;
 160        kesdtrack->GetXYZAt(0, xyz);
 161      }
 162      const Float_t kmass=kesdtrack->GetMass();
 163      const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
 164      const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
 165      const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
 166      const Float_t kp=TMath::Sqrt(kp2);
 167      const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30));
 168      const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
 169      //Double_t dx = xyz[0]-vertexpos[0];
 170      //Double_t dy = xyz[1]-vertexpos[1];
 171      //Float_t dca = TMath::Sqrt(dx*dx + dy*dy);
 172      //Float_t dz = xyz[2]-vertexpos[2];
 173      UInt_t index[6];
 174      const Int_t kncl=kesdtrack->GetITSclusters(index)
 175                       +kesdtrack->GetTPCclusters(NULL)
 176                       +kesdtrack->GetTRDclusters(NULL);
 177      if(IsAcceptedParticle(kpt,kphi,keta))
 178        fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kncl,kpt,kphi,keta);
 179
 180    } // loop over tracks
 181
 182   return kTRUE;
 183 }
 184
 185 /**********************************************************/
 186
 187 void AliJetParticlesReaderESD::Rewind()
 188 {
 189   //rewinds reading
 190
 191   if(fFile) delete fFile;
 192   if(fKeyIterator) delete fKeyIterator;
 193   fKeyIterator = 0;
 194   fFile = 0;
 195   fCurrentDir = 0;
 196   fNEventsRead = 0;
 197 }
 198
 199 /**********************************************************/
 200
 201 Int_t AliJetParticlesReaderESD::ReadNext()
 202 {
 203   //reads next event from fFile
 204
 205   do   // is OK even if 0 dirs specified,
 206     {  // in that case we try to read from "./"
 207       if (fFile == 0)
 208         {
 209           fFile = OpenFile(fCurrentDir);
 210           if (fFile == 0)
 211             {
 212               Error("ReadNext","Can't get fFile for dir no. %d",fCurrentDir);
 213               fCurrentDir++;
 214               continue;
 215             }
 216
 217           fCurrentEvent = 0;
 218           //fFile->GetListOfKeys()->Print();
 219
 220           if(fTree) delete fTree;
 221           fTree = dynamic_cast<TTree*>(fFile->Get("esdTree"));
 222           if(fTree)
 223             fTree->SetBranchAddress("ESD",&fESD);
 224           else
 225             fKeyIterator = new TIter(fFile->GetListOfKeys());
 226         }
 227
 228       if(fTree)
 229         {
 230           if(AliKalmanTrack::GetConvConst()<=0.)
 231             AliKalmanTrack::SetMagneticField(0.5);
 232           if(fCurrentEvent>=fTree->GetEntries())
 233             {
 234               fCurrentDir++;
 235               delete fTree;
 236               fTree = 0;
 237               delete fFile;
 238               fFile = 0;
 239               continue;
 240             }
 241           fTree->GetEvent(fCurrentEvent);
 242         }
 243       else
 244         { // "old" way via ESD objects stored in root file
 245           TKey* key = (TKey*)fKeyIterator->Next();
 246           if (key == 0)
 247             {
 248               fCurrentDir++;
 249               delete fKeyIterator;
 250               fKeyIterator = 0;
 251               delete fFile; //we have to assume there are no more ESD objects in the fFile
 252               fFile = 0;
 253               continue;
 254             }
 255           TString esdname = "ESD";
 256           esdname+=fCurrentEvent;
 257           if(fESD) delete fESD;
 258           fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
 259           if (fESD == 0)
 260             {
 261               Info("ReadNext","Can't find AliESD object named %s",esdname.Data());
 262               fCurrentDir++;
 263               delete fKeyIterator;
 264               fKeyIterator = 0;
 265               delete fFile;//we have to assume there is no more ESD objects in the fFile
 266               fFile = 0;
 267               continue;
 268             }
 269         }
 270       ReadESD(fESD);
 271       fCurrentEvent++;
 272       fNEventsRead++;
 273       return kTRUE;//success -> read one event
 274     }  while(fCurrentDir < GetNumberOfDirs());
 275       //end of loop over directories specified in fDirs Obj Array
 276   return kFALSE; //no more directories to read
 277 }
 278
 279 /**********************************************************/
 280
 281 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
 282 {
 283   //opens fFile with kine tree
 284
 285  const TString& dirname = GetDirName(n);
 286  if (dirname == "")
 287   {
 288    Error("OpenFiles","Can't get directory name");
 289    return 0;
 290   }
 291  TString filename = dirname +"/"+ fESDFileName;
 292  TFile *ret = TFile::Open(filename.Data());
 293  if (ret == 0)
 294   {
 295     Error("OpenFiles","Can't open fFile %s",filename.Data());
 296     return 0;
 297   }
 298  if (!ret->IsOpen())
 299   {
 300     Error("OpenFiles","Can't open fFile  %s",filename.Data());
 301     return 0;
 302   }
 303
 304  return ret;
 305 }
 306
 307 /**********************************************************/
 308
 309 void AliJetParticlesReaderESD::PrintESDtrack(const AliESDtrack *kesdtrack) const
 310 {
 311   ULong_t status=kesdtrack->GetStatus();
 312   cout << hex << status << dec << ": ";
 313   if((status & AliESDtrack::kITSin) == AliESDtrack::kITSin) cout << "ITSin ";
 314   if((status & AliESDtrack::kITSout) == AliESDtrack::kITSout) cout << "ITSout ";
 315   if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) cout << "ITSrefit ";
 316   if((status & AliESDtrack::kITSpid) == AliESDtrack::kITSpid) cout << "ITSpid ";
 317
 318   if((status & AliESDtrack::kTPCin)  == AliESDtrack::kTPCin) cout << "TPCin ";
 319   if((status & AliESDtrack::kTPCout) == AliESDtrack::kTPCout) cout << "TPCout ";
 320   if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) cout << "TPCrefit ";
 321   if((status & AliESDtrack::kTPCpid) == AliESDtrack::kTPCpid) cout << "TPCpid ";
 322
 323   if((status & AliESDtrack::kTRDin) == AliESDtrack::kTRDin) cout << "TRDin ";
 324   if((status & AliESDtrack::kTRDout) == AliESDtrack::kTRDout) cout << "TRDout ";
 325   if((status & AliESDtrack::kTRDrefit) == AliESDtrack::kTRDrefit) cout << "TRDrefit ";
 326   if((status & AliESDtrack::kTRDpid) == AliESDtrack::kTRDpid) cout << "TRDpid ";
 327   if((status & AliESDtrack::kTRDbackup) == AliESDtrack::kTRDbackup) cout << "TRDbackup ";
 328   if((status & AliESDtrack::kTRDStop) == AliESDtrack::kTRDStop) cout << "TRDstop ";
 329
 330   if((status & AliESDtrack::kTOFin) == AliESDtrack::kTOFin) cout << "TOFin ";
 331   if((status & AliESDtrack::kTOFout) == AliESDtrack::kTOFout) cout << "TOFout ";
 332   if((status & AliESDtrack::kTOFrefit) == AliESDtrack::kTOFrefit) cout << "TOFrefit ";
 333   if((status & AliESDtrack::kTOFpid) == AliESDtrack::kTOFpid) cout << "TOFpid ";
 334
 335   if((status & AliESDtrack::kPHOSpid) == AliESDtrack::kPHOSpid) cout << "PHOSpid ";
 336   if((status & AliESDtrack::kRICHpid) == AliESDtrack::kRICHpid) cout << "RICHpid ";
 337   if((status & AliESDtrack::kEMCALpid) == AliESDtrack::kEMCALpid) cout << "EMCALpid ";
 338   if((status & AliESDtrack::kESDpid) == AliESDtrack::kESDpid) cout << "ESDpid ";
 339   if((status & AliESDtrack::kTIME) == AliESDtrack::kTIME) cout << "TIME ";
 340   cout << endl;
 341   /*
 342   cout << kesdtrack->GetConstrainedChi2()
 343        << " " << kesdtrack->GetITSchi2()
 344        << " " << kesdtrack->GetTPCchi2()
 345        << " " << kesdtrack->GetTRDchi2()<< endl;
 346   */
 347 }