HBTAN/AliHBTReaderESD.cxx

   1 #include "AliHBTReaderESD.h"
   2 //____________________________________________________________________
   3 //////////////////////////////////////////////////////////////////////
   4 //                                                                  //
   5 // class AliHBTReaderESD                                            //
   6 //                                                                  //
   7 // reader for ALICE Event Summary Data (ESD).                       //
   8 //                                                                  //
   9 // Piotr.Skowronski@cern.ch                                         //
  10 //                                                                  //
  11 //////////////////////////////////////////////////////////////////////
  12
  13 #include <TPDGCode.h>
  14 #include <TString.h>
  15 #include <TObjString.h>
  16 #include <TTree.h>
  17 #include <TFile.h>
  18 #include <TKey.h>
  19 #include <TParticle.h>
  20 #include <TH1.h>
  21
  22 #include <AliRun.h>
  23 #include <AliRunLoader.h>
  24 #include <AliStack.h>
  25 #include <AliESDtrack.h>
  26 #include <AliESD.h>
  27
  28 #include "AliHBTRun.h"
  29 #include "AliHBTEvent.h"
  30 #include "AliHBTParticle.h"
  31 #include "AliHBTParticleCut.h"
  32 #include "AliHBTTrackPoints.h"
  33 #include "AliHBTClusterMap.h"
  34
  35 ClassImp(AliHBTReaderESD)
  36
  37 AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
  38  fESDFileName(esdfilename),
  39  fGAlFileName(galfilename),
  40  fFile(0x0),
  41  fRunLoader(0x0),
  42  fKeyIterator(0x0),
  43  fReadParticles(kFALSE),
  44  fNTrackPoints(0),
  45  fdR(0.0),
  46  fClusterMap(kFALSE),
  47  fNTPCClustMin(0),
  48  fNTPCClustMax(150),
  49  fTPCChi2PerClustMin(0.0),
  50  fTPCChi2PerClustMax(10e5),
  51  fC00Min(0.0),
  52  fC00Max(10e5),
  53  fC11Min(0.0),
  54  fC11Max(10e5),
  55  fC22Min(0.0),
  56  fC22Max(10e5),
  57  fC33Min(0.0),
  58  fC33Max(10e5),
  59  fC44Min(0.0),
  60  fC44Max(10e5),
  61  fTPCC00Min(0.0),
  62  fTPCC00Max(10e5),
  63  fTPCC11Min(0.0),
  64  fTPCC11Max(10e5),
  65  fTPCC22Min(0.0),
  66  fTPCC22Max(10e5),
  67  fTPCC33Min(0.0),
  68  fTPCC33Max(10e5),
  69  fTPCC44Min(0.0),
  70  fTPCC44Max(10e5)
  71
  72 {
  73   //cosntructor
  74   if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
  75     Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
  76 }
  77 /********************************************************************/
  78
  79 AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
  80  AliHBTReader(dirs),
  81  fESDFileName(esdfilename),
  82  fGAlFileName(galfilename),
  83  fFile(0x0),
  84  fRunLoader(0x0),
  85  fKeyIterator(0x0),
  86  fReadParticles(kFALSE),
  87  fNTrackPoints(0),
  88  fdR(0.0),
  89  fClusterMap(kFALSE),
  90  fNTPCClustMin(0),
  91  fNTPCClustMax(150),
  92  fTPCChi2PerClustMin(0.0),
  93  fTPCChi2PerClustMax(10e5),
  94  fC00Min(0.0),
  95  fC00Max(10e5),
  96  fC11Min(0.0),
  97  fC11Max(10e5),
  98  fC22Min(0.0),
  99  fC22Max(10e5),
 100  fC33Min(0.0),
 101  fC33Max(10e5),
 102  fC44Min(0.0),
 103  fC44Max(10e5),
 104  fTPCC00Min(0.0),
 105  fTPCC00Max(10e5),
 106  fTPCC11Min(0.0),
 107  fTPCC11Max(10e5),
 108  fTPCC22Min(0.0),
 109  fTPCC22Max(10e5),
 110  fTPCC33Min(0.0),
 111  fTPCC33Max(10e5),
 112  fTPCC44Min(0.0),
 113  fTPCC44Max(10e5)
 114 {
 115   //cosntructor
 116   if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
 117     Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
 118 }
 119 /********************************************************************/
 120
 121 AliHBTReaderESD::~AliHBTReaderESD()
 122 {
 123  //desctructor
 124   delete fRunLoader;
 125   delete fKeyIterator;
 126   delete fFile;
 127 }
 128 /**********************************************************/
 129 Int_t AliHBTReaderESD::ReadNext()
 130 {
 131 //reads next event from fFile
 132   //fRunLoader is for reading Kine
 133
 134   if (AliHBTParticle::GetDebug())
 135     Info("ReadNext","Entered");
 136
 137   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
 138   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
 139
 140   fParticlesEvent->Reset();
 141   fTracksEvent->Reset();
 142
 143   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
 144    {
 145      if (fFile == 0x0)
 146       {
 147        fFile = OpenFile(fCurrentDir);//rl is opened here
 148        if (fFile == 0x0)
 149          {
 150            Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
 151            fCurrentDir++;
 152            continue;
 153          }
 154        fCurrentEvent = 0;
 155        fKeyIterator = new TIter(fFile->GetListOfKeys());
 156 //       fFile->Dump();
 157 //       fFile->GetListOfKeys()->Print();
 158       }
 159      TKey* key = (TKey*)fKeyIterator->Next();
 160      if (key == 0x0)
 161       {
 162         if (AliHBTParticle::GetDebug() > 2 )
 163           {
 164             Info("ReadNext","No more keys.");
 165           }
 166         fCurrentDir++;
 167         delete fKeyIterator;
 168         fKeyIterator = 0x0;
 169         delete fFile;//we have to assume there is no more ESD objects in the fFile
 170         fFile = 0x0;
 171         delete fRunLoader;
 172         fRunLoader = 0x0;
 173         continue;
 174       }
 175      //try to read
 176
 177
 178 //     TObject* esdobj = key->ReadObj();
 179 //     if (esdobj == 0x0)
 180 //      {
 181 //        if (AliHBTParticle::GetDebug() > 2 )
 182 //          {
 183 //            Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
 184 //            key->Dump();
 185 //          }
 186 //        continue;
 187 //      }
 188 //     esdobj->Dump();
 189 //     AliESD* esd = dynamic_cast<AliESD*>(esdobj);
 190
 191      TString esdname = "ESD";
 192      esdname+=fCurrentEvent;
 193      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
 194      if (esd == 0x0)
 195       {
 196 //        if (AliHBTParticle::GetDebug() > 2 )
 197 //          {
 198 //            Info("ReadNext","This key is not an AliESD object %s",key->GetName());
 199 //          }
 200         if (AliHBTParticle::GetDebug() > 2 )
 201           {
 202             Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
 203           }
 204         fCurrentDir++;
 205         delete fKeyIterator;
 206         fKeyIterator = 0x0;
 207         delete fFile;//we have to assume there is no more ESD objects in the fFile
 208         fFile = 0x0;
 209         delete fRunLoader;
 210         fRunLoader = 0x0;
 211         continue;
 212       }
 213
 214      ReadESD(esd);
 215
 216      fCurrentEvent++;
 217      fNEventsRead++;
 218      delete esd;
 219      return 0;//success -> read one event
 220    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
 221
 222   return 1; //no more directories to read
 223 }
 224 /**********************************************************/
 225
 226 Int_t AliHBTReaderESD::ReadESD(AliESD* esd)
 227 {
 228   //****** Tentative particle type "concentrations"
 229   static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
 230
 231   Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
 232   Double_t w[kNSpecies];
 233   Double_t mom[3];//momentum
 234   Double_t pos[3];//position
 235   Double_t vertexpos[3];//vertex position
 236   //Reads one ESD
 237   if (esd == 0x0)
 238    {
 239      Error("ReadESD","ESD is NULL");
 240      return 1;
 241    }
 242
 243   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
 244   if (pdgdb == 0x0)
 245    {
 246      Error("ReadESD","Can not get PDG Database Instance.");
 247      return 1;
 248    }
 249
 250   Float_t mf = esd->GetMagneticField();
 251
 252   if ( (mf == 0.0) && (fNTrackPoints > 0) )
 253    {
 254       Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
 255       return 1;
 256    }
 257
 258   AliStack* stack = 0x0;
 259   if (fReadParticles && fRunLoader)
 260    {
 261      fRunLoader->GetEvent(fCurrentEvent);
 262      stack = fRunLoader->Stack();
 263    }
 264
 265   const AliESDVertex* vertex = esd->GetVertex();
 266   if (vertex == 0x0)
 267    {
 268      Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
 269      vertexpos[0] = 0.0;
 270      vertexpos[1] = 0.0;
 271      vertexpos[2] = 0.0;
 272    }
 273   else
 274    {
 275      vertex->GetXYZ(vertexpos);
 276    }
 277
 278   if (AliHBTParticle::GetDebug() > 0)
 279    {
 280      Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]);
 281    }
 282
 283   Info("ReadESD","Reading Event %d",fCurrentEvent);
 284
 285   Int_t ntr = esd->GetNumberOfTracks();
 286   Info("ReadESD","Found %d tracks.",ntr);
 287   for (Int_t i = 0;i<ntr; i++)
 288    {
 289      AliESDtrack *esdtrack = esd->GetTrack(i);
 290      if (esdtrack == 0x0)
 291       {
 292         Error("Next","Can not get track %d", i);
 293         continue;
 294       }
 295
 296      //if (esdtrack->HasVertexParameters() == kFALSE)
 297      if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
 298       {
 299         if (AliHBTParticle::GetDebug() > 2)
 300           Info("ReadNext","Particle skipped: Data at vertex not available.");
 301         continue;
 302       }
 303
 304      if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
 305       {
 306         if (AliHBTParticle::GetDebug() > 2)
 307           Info("ReadNext","Particle skipped: PID BIT is not set.");
 308         continue;
 309       }
 310
 311
 312      Double_t extx;
 313      Double_t extp[5];
 314      esdtrack->GetConstrainedExternalParameters(extx,extp);
 315      if (extp[4] == 0.0)
 316       {
 317         if (AliHBTParticle::GetDebug() > 2)
 318           Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
 319         continue;
 320       }
 321      esdtrack->GetESDpid(pidtable);
 322      //esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]);
 323      esdtrack->GetConstrainedPxPyPz(mom);
 324      //esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]);
 325      esdtrack->GetConstrainedXYZ(pos);
 326      pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
 327      pos[1] -= vertexpos[1];
 328      pos[2] -= vertexpos[2];
 329
 330      Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
 331
 332      //Particle from kinematics
 333      AliHBTParticle* particle = 0;
 334      Bool_t keeppart = kFALSE;
 335      if ( fReadParticles && stack  )
 336       {
 337         if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track
 338         TParticle *p = stack->Particle(esdtrack->GetLabel());
 339         if (p==0x0)
 340          {
 341            Error("ReadNext","Can not find track with such label.");
 342            continue;
 343          }
 344 //           if(p->GetPdgCode()<0) charge = -1;
 345         particle = new AliHBTParticle(*p,i);
 346
 347       }
 348
 349      //Here we apply Bayes' formula
 350      Double_t rc=0.;
 351      for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
 352      if (rc==0.0)
 353       {
 354         if (AliHBTParticle::GetDebug() > 2)
 355           Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
 356         continue;
 357       }
 358
 359      for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
 360
 361      if (AliHBTParticle::GetDebug() > 4)
 362       {
 363         Info("ReadNext","###########################################################################");
 364         Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
 365         Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
 366         TString msg("Pid list got from track:");
 367         for (Int_t s = 0;s<kNSpecies;s++)
 368          {
 369            msg+="\n    ";
 370            msg+=s;
 371            msg+="(";
 372            msg+=charge*GetSpeciesPdgCode((ESpecies)s);
 373            msg+="): ";
 374            msg+=w[s];
 375            msg+=" (";
 376            msg+=pidtable[s];
 377            msg+=")";
 378          }
 379         Info("ReadNext","%s",msg.Data());
 380       }//if (AliHBTParticle::GetDebug()>4)
 381
 382       AliHBTTrackPoints* tpts = 0x0;
 383       if (fNTrackPoints > 0)
 384        {
 385          tpts = new AliHBTTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
 386        }
 387
 388       AliHBTClusterMap* cmap = 0x0;
 389       if (  fClusterMap )
 390        {
 391          cmap = new AliHBTClusterMap(esdtrack);
 392        }
 393
 394      for (Int_t s = 0; s<kNSpecies; s++)
 395       {
 396         Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
 397         Float_t pp = w[s];
 398         if (pp == 0.0)
 399          {
 400            if ( AliHBTParticle::GetDebug() > 5 )
 401              Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
 402            continue;
 403          }
 404
 405         if(Pass(pdgcode))
 406          {
 407            if ( AliHBTParticle::GetDebug() > 5 )
 408              Info("ReadNext","PID (%d) did not pass the cut.",pdgcode);
 409            continue; //check if we are intersted with particles of this type
 410          }
 411
 412         Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
 413         Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
 414
 415         AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
 416                                                    mom[0], mom[1], mom[2], tEtot,
 417                                                    pos[0], pos[1], pos[2], 0.);
 418         //copy probabilitis of other species (if not zero)
 419         for (Int_t k = 0; k<kNSpecies; k++)
 420          {
 421            if (k == s) continue;
 422            if (w[k] == 0.0) continue;
 423            track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
 424          }
 425
 426         if(Pass(track))//check if meets all criteria of any of our cuts
 427                        //if it does not delete it and take next good track
 428          {
 429            if ( AliHBTParticle::GetDebug() > 4 )
 430              Info("ReadNext","Track did not pass the cut");
 431            delete track;
 432            continue;
 433          }
 434
 435         //Single Particle cuts on cluster map and track points rather do not have sense
 436         if (tpts)
 437          {
 438            track->SetTrackPoints(tpts);
 439          }
 440
 441         if (cmap)
 442          {
 443            track->SetClusterMap(cmap);
 444          }
 445
 446         fTracksEvent->AddParticle(track);
 447         if (particle) fParticlesEvent->AddParticle(particle);
 448         keeppart = kTRUE;
 449
 450         if (AliHBTParticle::GetDebug() > 4 )
 451          {
 452            Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
 453            track->Print();
 454            if (particle) particle->Print();
 455            Info("ReadNext","\n----------------------------------------------\n");
 456          }
 457       }//for (Int_t s = 0; s<kNSpecies; s++)
 458
 459      if (keeppart == kFALSE)
 460       {
 461         delete particle;//particle was not stored in event
 462         delete tpts;
 463         delete cmap;
 464       }
 465
 466    }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
 467
 468   Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
 469          fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
 470          fNEventsRead,fCurrentEvent,fCurrentDir);
 471   fTrackCounter->Fill(fTracksEvent->GetNumberOfParticles());
 472   return 0;
 473 }
 474
 475 /**********************************************************/
 476
 477 void AliHBTReaderESD::Rewind()
 478 {
 479   //rewinds reading
 480   delete fKeyIterator;
 481   delete fFile;
 482   fFile = 0x0;
 483   fKeyIterator = 0x0;
 484   delete fRunLoader;
 485   fRunLoader = 0x0;
 486   fCurrentDir = 0;
 487   fNEventsRead = 0;
 488   if (fTrackCounter) fTrackCounter->Reset();
 489 }
 490 /**********************************************************/
 491
 492 TFile* AliHBTReaderESD::OpenFile(Int_t n)
 493 {
 494 //opens fFile with kine tree
 495
 496  const TString& dirname = GetDirName(n);
 497  if (dirname == "")
 498   {
 499    Error("OpenFiles","Can not get directory name");
 500    return 0x0;
 501   }
 502  TString filename = dirname +"/"+ fESDFileName;
 503  TFile *ret = TFile::Open(filename.Data());
 504
 505  if ( ret == 0x0)
 506   {
 507     Error("OpenFiles","Can't open fFile %s",filename.Data());
 508     return 0x0;
 509   }
 510  if (!ret->IsOpen())
 511   {
 512     Error("OpenFiles","Can't open fFile  %s",filename.Data());
 513     return 0x0;
 514   }
 515
 516  if (fReadParticles )
 517   {
 518    fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
 519    if (fRunLoader == 0x0)
 520     {
 521       Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
 522       delete ret;
 523       return 0x0;
 524     }
 525
 526    fRunLoader->LoadHeader();
 527    if (fRunLoader->LoadKinematics())
 528     {
 529       Error("Next","Error occured while loading kinematics.");
 530       delete fRunLoader;
 531       delete ret;
 532       return 0x0;
 533     }
 534   }
 535
 536  return ret;
 537 }
 538 /**********************************************************/
 539
 540 Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
 541 {
 542   //returns pdg code from the PID index
 543   //ask jura about charge
 544   switch (spec)
 545    {
 546      case kESDElectron:
 547        return kPositron;
 548        break;
 549      case kESDMuon:
 550        return kMuonPlus;
 551        break;
 552      case kESDPion:
 553        return kPiPlus;
 554        break;
 555      case kESDKaon:
 556        return kKPlus;
 557        break;
 558      case kESDProton:
 559        return kProton;
 560        break;
 561      default:
 562        ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
 563        break;
 564    }
 565   return 0;
 566 }
 567 /********************************************************************/
 568
 569 void AliHBTReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
 570 {
 571  //sets range of Number Of Clusters that tracks have to have
 572  fNTPCClustMin = min;
 573  fNTPCClustMax = max;
 574 }
 575 /********************************************************************/
 576
 577 void AliHBTReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
 578 {
 579   //sets range of Chi2 per Cluster
 580   fTPCChi2PerClustMin = min;
 581   fTPCChi2PerClustMax = max;
 582 }
 583 /********************************************************************/
 584
 585 void AliHBTReaderESD::SetC00Range(Float_t min, Float_t max)
 586 {
 587  //Sets range of C00 parameter of covariance matrix of the track
 588  //it defines uncertainty of the momentum
 589   fC00Min = min;
 590   fC00Max = max;
 591 }
 592 /********************************************************************/
 593
 594 void AliHBTReaderESD::SetC11Range(Float_t min, Float_t max)
 595 {
 596  //Sets range of C11 parameter of covariance matrix of the track
 597  //it defines uncertainty of the momentum
 598   fC11Min = min;
 599   fC11Max = max;
 600 }
 601 /********************************************************************/
 602
 603 void AliHBTReaderESD::SetC22Range(Float_t min, Float_t max)
 604 {
 605  //Sets range of C22 parameter of covariance matrix of the track
 606  //it defines uncertainty of the momentum
 607   fC22Min = min;
 608   fC22Max = max;
 609 }
 610 /********************************************************************/
 611
 612 void AliHBTReaderESD::SetC33Range(Float_t min, Float_t max)
 613 {
 614  //Sets range of C33 parameter of covariance matrix of the track
 615  //it defines uncertainty of the momentum
 616   fC33Min = min;
 617   fC33Max = max;
 618 }
 619 /********************************************************************/
 620
 621 void AliHBTReaderESD::SetC44Range(Float_t min, Float_t max)
 622 {
 623  //Sets range of C44 parameter of covariance matrix of the track
 624  //it defines uncertainty of the momentum
 625   fC44Min = min;
 626   fC44Max = max;
 627 }