PWG2/RESONANCES/AliRsnAnalysisMonitorPairTask.cxx

   1 //
   2 // Implementation file for implementation of data analysis aft 900 GeV
   3 //
   4 // Author: A. Pulvirenti
   5 //
   6
   7 #include "Riostream.h"
   8 #include <iomanip>
   9
  10 #include "TH1.h"
  11 #include "TTree.h"
  12 #include "TParticle.h"
  13 #include "TRandom.h"
  14 #include "TLorentzVector.h"
  15
  16 #include "AliLog.h"
  17 #include "AliESDpid.h"
  18 #include "AliESDEvent.h"
  19 #include "AliESDVertex.h"
  20 #include "AliESDtrack.h"
  21 #include "AliStack.h"
  22 #include "AliMCEvent.h"
  23 #include "AliTOFT0maker.h"
  24 #include "AliTOFcalib.h"
  25 #include "AliCDBManager.h"
  26 #include "AliITSPIDResponse.h"
  27 #include "AliRsnMonitorTrack.h"
  28 #include "AliRsnDaughter.h"
  29 #include "AliRsnEvent.h"
  30
  31 #include "AliRsnAnalysisMonitorPairTask.h"
  32
  33 ClassImp(AliRsnAnalysisMonitorPairTask)
  34
  35 //__________________________________________________________________________________________________
  36 AliRsnAnalysisMonitorPairTask::AliRsnAnalysisMonitorPairTask(const char *name) :
  37   AliAnalysisTaskSE(name),
  38   fOut(0x0),
  39   fInvMass(0.0),
  40   fRangeMin(0.0),
  41   fRangeMax(100.0),
  42   fESDpid(0x0),
  43   fTOFmaker(0x0),
  44   fTOFcalib(0x0),
  45   fTOFcalibrateESD(kFALSE),
  46   fTOFcorrectTExp(kFALSE),
  47   fTOFuseT0(kFALSE),
  48   fTOFtuneMC(kFALSE),
  49   fTOFresolution(0.0),
  50   fEventCuts("eventCuts", AliRsnCut::kEvent),
  51   fTrackCuts("trackCuts", AliRsnCut::kDaughter)
  52
  53 {
  54 //
  55 // Constructor
  56 //
  57
  58   DefineOutput(1, TTree::Class());
  59 }
  60
  61 //__________________________________________________________________________________________________
  62 AliRsnAnalysisMonitorPairTask::AliRsnAnalysisMonitorPairTask(const AliRsnAnalysisMonitorPairTask& copy) :
  63   AliAnalysisTaskSE(copy),
  64   fOut(0x0),
  65   fInvMass(0.0),
  66   fRangeMin(copy.fRangeMin),
  67   fRangeMax(copy.fRangeMax),
  68   fESDpid(0x0),
  69   fTOFmaker(0x0),
  70   fTOFcalib(0x0),
  71   fTOFcalibrateESD(kFALSE),
  72   fTOFcorrectTExp(kFALSE),
  73   fTOFuseT0(kFALSE),
  74   fTOFtuneMC(kFALSE),
  75   fTOFresolution(0.0),
  76   fEventCuts(copy.fEventCuts),
  77   fTrackCuts(copy.fTrackCuts)
  78 {
  79 //
  80 // Copy constructor
  81 //
  82 }
  83
  84 //__________________________________________________________________________________________________
  85 AliRsnAnalysisMonitorPairTask& AliRsnAnalysisMonitorPairTask::operator=(const AliRsnAnalysisMonitorPairTask& copy)
  86 {
  87 //
  88 // Assignment operator
  89 //
  90
  91   fMass[0] = copy.fMass[0];
  92   fMass[1] = copy.fMass[1];
  93
  94   fRangeMin = copy.fRangeMin;
  95   fRangeMax = copy.fRangeMax;
  96
  97   fTOFcalibrateESD = copy.fTOFcalibrateESD;
  98   fTOFcorrectTExp = copy.fTOFcorrectTExp;
  99   fTOFuseT0 = copy.fTOFuseT0;
 100   fTOFtuneMC = copy.fTOFtuneMC;
 101   fTOFresolution = copy.fTOFresolution;
 102
 103   return (*this);
 104 }
 105
 106 //__________________________________________________________________________________________________
 107 AliRsnAnalysisMonitorPairTask::~AliRsnAnalysisMonitorPairTask()
 108 {
 109 //
 110 // Destructor
 111 //
 112
 113   if (fOut)    delete fOut;
 114   if (fESDpid) delete fESDpid;
 115   if (fTrack)  delete fTrack;
 116 }
 117
 118 //__________________________________________________________________________________________________
 119 void AliRsnAnalysisMonitorPairTask::UserCreateOutputObjects()
 120 {
 121 //
 122 // Create the output data container
 123 //
 124
 125   // setup TPC response
 126   fESDpid = new AliESDpid;
 127   fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0], fTPCpar[1], fTPCpar[2], fTPCpar[3], fTPCpar[4]);
 128
 129   // setup TOF maker & calibration
 130   fTOFcalib = new AliTOFcalib;
 131   fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
 132   fTOFmaker->SetTimeResolution(fTOFresolution);
 133
 134   // create output branch object
 135   fTrack[0] = new AliRsnMonitorTrack;
 136   fTrack[1] = new AliRsnMonitorTrack;
 137
 138   // create output tree
 139   OpenFile(1);
 140   fOut = new TTree("rsnPairMonitor", "Informations on pairs for cut checking");
 141   fOut->Branch("track1", "AliRsnMonitorTrack", &fTrack[0]);
 142   fOut->Branch("track2", "AliRsnMonitorTrack", &fTrack[1]);
 143   fOut->Branch("minv"  , &fInvMass, "minv/F");
 144 }
 145
 146 //__________________________________________________________________________________________________
 147 void AliRsnAnalysisMonitorPairTask::UserExec(Option_t *)
 148 {
 149 //
 150 // Main execution function.
 151 // Fills the fHEvents data member with the following legenda:
 152 // 0 -- event OK, prim vertex with tracks
 153 // 1 -- event OK, prim vertex with SPD
 154 // 2 -- event OK but vz large
 155 // 3 -- event bad
 156 //
 157
 158   // retrieve ESD event and related stack (if available)
 159   AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
 160   AliStack    *stack = 0x0;
 161
 162   // skip NULL events
 163   if (!esd) return;
 164   if (fMCEvent) stack = fMCEvent->Stack();
 165
 166   // create interface objects to AliRsnEvent to check event cuts
 167   AliRsnEvent event;
 168   event.SetRef(esd);
 169   event.SetRefMC(fMCEvent);
 170   if (!fEventCuts.IsSelected(&event)) return;
 171
 172   // check the event
 173   Int_t type = EventEval(esd);
 174
 175   // if processable, then process it
 176   if      (type == 0) ProcessESD(esd, esd->GetPrimaryVertexTracks());
 177   else if (type == 1) ProcessESD(esd, esd->GetPrimaryVertexSPD());
 178   else return;
 179
 180   // update histogram container
 181   PostData(1, fOut);
 182 }
 183
 184 //__________________________________________________________________________________________________
 185 void AliRsnAnalysisMonitorPairTask::Terminate(Option_t *)
 186 {
 187 //
 188 // Terminate
 189 //
 190 }
 191
 192 //__________________________________________________________________________________________________
 193 Int_t AliRsnAnalysisMonitorPairTask::EventEval(AliESDEvent *esd)
 194 {
 195 //
 196 // Checks if the event is good for analysis.
 197 // Returns:
 198 // ---> 0 if a good primary vertex with tracks was found,
 199 // ---> 1 if a good SPD primary vertex was found
 200 // ---> 2 otherwise (event to be rejected)
 201 // In any case, adds an entry to the TTree, to keep trace of all events.
 202 //
 203
 204   // get the best primary vertex:
 205   // first try that with tracks, then the SPD one
 206   const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
 207   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
 208   if(vTrk->GetNContributors() > 0)
 209   {
 210     return 0;
 211   }
 212   else if (vSPD->GetNContributors() > 0)
 213   {
 214     return 1;
 215   }
 216   else
 217   {
 218     return 2;
 219   }
 220 }
 221
 222 //__________________________________________________________________________________________________
 223 void AliRsnAnalysisMonitorPairTask::ProcessESD(AliESDEvent *esd, const AliESDVertex *v)
 224 {
 225 //
 226 // Process the ESD container, to read all tracks and copy their useful values.
 227 // All info are stored into an AliRsnMonitorTrack object and saved into the
 228 // TClonesArray which is one of the branches of the output TTree.
 229 //
 230
 231   // TOF stuff #1: init OCDB
 232   Int_t run = esd->GetRunNumber();
 233   AliCDBManager *cdb = AliCDBManager::Instance();
 234   cdb->SetDefaultStorage("raw://");
 235   cdb->SetRun(run);
 236   // TOF stuff #2: init calibration
 237   fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 238   fTOFcalib->Init();
 239   // TOF stuff #3: calibrate
 240   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
 241   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
 242   if (fTOFuseT0)
 243   {
 244     fTOFmaker->ComputeT0TOF(esd);
 245     fTOFmaker->ApplyT0TOF(esd);
 246     fESDpid->MakePID(esd, kFALSE, 0.);
 247   }
 248
 249   // loop on all tracks
 250   Int_t          i0 , i1, ch0, ch1, ntracks = esd->GetNumberOfTracks();
 251   TLorentzVector v0, v1, sum;
 252
 253   // two nested loops on tracks
 254   for (i0 = 0; i0 < ntracks; i0++)
 255   {
 256     if (!ProcessTrack(0, i0, esd, v)) continue;
 257     ch0 = fTrack[0]->Charge();
 258     v0.SetXYZM(fTrack[0]->PrecX(), fTrack[0]->PrecY(), fTrack[0]->PrecZ(), fMass[0]);
 259
 260     for (i1 = i0 + 1; i1 < ntracks; i1++)
 261     {
 262       if (!ProcessTrack(1, i1, esd, v)) continue;
 263       ch1 = fTrack[1]->Charge();
 264
 265       // skip like-sign pairs
 266       if (ch1 == ch0) continue;
 267
 268       // check invmass range
 269       v1.SetXYZM(fTrack[1]->PrecX(), fTrack[1]->PrecY(), fTrack[1]->PrecZ(), fMass[1]);
 270       sum = v0 + v1;
 271       fInvMass = sum.M();
 272       if (fInvMass < fRangeMin || fInvMass > fRangeMax) continue;
 273
 274       // if here, add an entry to the tree
 275       fOut->Fill();
 276     }
 277   }
 278 }
 279
 280 //__________________________________________________________________________________________________
 281 Bool_t AliRsnAnalysisMonitorPairTask::ProcessTrack(Int_t myIndex, Int_t esdIndex, AliESDEvent *esd, const AliESDVertex *v)
 282 {
 283 //
 284 // Process a single track, and stores all its info into one of the two branch objects
 285 // defined by the first index, while the second chooses the track in the owner ESD
 286 //
 287
 288   if (myIndex < 0 || myIndex > 1) return kFALSE;
 289
 290   // ITS stuff #1
 291   // create the response function and initialize it to MC or not
 292   // depending if the AliStack object is there or not
 293   AliStack *stack = 0x0;
 294   if (fMCEvent) stack = fMCEvent->Stack();
 295   Bool_t isMC = (stack != 0x0);
 296   AliITSPIDResponse itsrsp(isMC);
 297
 298   // create interfacr objects
 299   AliRsnEvent    event;
 300   event.SetRef(esd);
 301   event.SetRefMC(fMCEvent);
 302   AliRsnDaughter daughter;
 303   AliESDtrack *track = esd->GetTrack(esdIndex);
 304   event.SetDaughter(daughter, esdIndex, AliRsnDaughter::kTrack);
 305
 306   // skip NULL pointers, kink daughters and tracks which
 307   // cannot be propagated to primary vertex
 308   if (!track) return kFALSE;
 309   if ((Int_t)track->GetKinkIndex(0) > 0) return kFALSE;
 310   if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) return kFALSE;
 311
 312   // useful variables
 313   Int_t     k, nITS;
 314   Bool_t    isTPC, isITSSA, isTOF;
 315   Float_t   b[2], bCov[3];
 316   Double_t  time[10];
 317
 318   // reset the output object
 319   // 'usable' flag will need to be set to 'ok'
 320   fTrack[myIndex]->Reset();
 321
 322   // check cuts
 323   fTrack[myIndex]->CutsPassed() = fTrackCuts.IsSelected(&daughter);
 324
 325   // get MC info if possible
 326   if (stack) fTrack[myIndex]->AdoptMC(TMath::Abs(track->GetLabel()), stack);
 327
 328   // copy general info
 329   fTrack[myIndex]->Status() = (UInt_t)track->GetStatus();
 330   fTrack[myIndex]->Length() = (Double_t)track->GetIntegratedLength();
 331   fTrack[myIndex]->Charge() = (Int_t)track->Charge();
 332   fTrack[myIndex]->PrecX()  = (Double_t)track->Px();
 333   fTrack[myIndex]->PrecY()  = (Double_t)track->Py();
 334   fTrack[myIndex]->PrecZ()  = (Double_t)track->Pz();
 335
 336   // evaluate some flags from the status to decide what to do next in some points
 337   isTPC   = ((fTrack[myIndex]->Status() & AliESDtrack::kTPCin)  != 0);
 338   isITSSA = ((fTrack[myIndex]->Status() & AliESDtrack::kTPCin)  == 0 && (fTrack[myIndex]->Status() & AliESDtrack::kITSrefit) != 0 && (fTrack[myIndex]->Status() & AliESDtrack::kITSpureSA) == 0 && (fTrack[myIndex]->Status() & AliESDtrack::kITSpid) != 0);
 339   isTOF   = ((fTrack[myIndex]->Status() & AliESDtrack::kTOFout) != 0 && (fTrack[myIndex]->Status() & AliESDtrack::kTIME) != 0);
 340
 341   // accept only tracks which are TPC+ITS or ITS standalone
 342   if (isITSSA)
 343   {
 344     fTrack[myIndex]->ITSsa() = kTRUE;
 345     fTrack[myIndex]->TOFok() = kFALSE;
 346   }
 347   else if (isTPC)
 348   {
 349     fTrack[myIndex]->ITSsa() = kFALSE;
 350     fTrack[myIndex]->TOFok() = isTOF;
 351   }
 352   else
 353     return kFALSE;
 354
 355   // get DCA to primary vertex
 356   track->GetImpactParameters(b, bCov);
 357   fTrack[myIndex]->DCAr() = (Double_t)b[0];
 358   fTrack[myIndex]->DCAz() = (Double_t)b[1];
 359
 360   // get ITS info
 361   for (k = 0; k < 6; k++)
 362   {
 363     fTrack[myIndex]->ITSmap(k) = track->HasPointOnITSLayer(k);
 364   }
 365   if (isITSSA)
 366   {
 367     fTrack[myIndex]->ITSchi2() = track->GetITSchi2();
 368     fTrack[myIndex]->ITSsignal() = track->GetITSsignal();
 369     nITS = fTrack[myIndex]->SSDcount() + fTrack[myIndex]->SDDcount();
 370     for (k = 0; k < AliPID::kSPECIES; k++)
 371     {
 372       fTrack[myIndex]->ITSnsigma(k) = itsrsp.GetNumberOfSigmas(fTrack[myIndex]->Prec(), fTrack[myIndex]->ITSsignal(), (AliPID::EParticleType)k, nITS, kTRUE);
 373     }
 374   }
 375
 376   // get TPC info
 377   if (isTPC)
 378   {
 379     fTrack[myIndex]->TPCcount()  = (Int_t)track->GetTPCclusters(0);
 380     fTrack[myIndex]->TPCchi2()   = (Double_t)track->GetTPCchi2();
 381     fTrack[myIndex]->TPCsignal() = (Double_t)track->GetTPCsignal();
 382     fTrack[myIndex]->PtpcX()     = fTrack[myIndex]->PtpcY() = fTrack[myIndex]->PtpcZ() = 1E10;
 383     if (track->GetInnerParam())
 384     {
 385       fTrack[myIndex]->PtpcX() = track->GetInnerParam()->Px();
 386       fTrack[myIndex]->PtpcY() = track->GetInnerParam()->Py();
 387       fTrack[myIndex]->PtpcZ() = track->GetInnerParam()->Pz();
 388       for (k = 0; k < AliPID::kSPECIES; k++)
 389       {
 390         fTrack[myIndex]->TPCnsigma(k) = fESDpid->NumberOfSigmasTPC(track, (AliPID::EParticleType)k);
 391       }
 392     }
 393   }
 394
 395   // get TOF info
 396   if (isTOF)
 397   {
 398     track->GetIntegratedTimes(time);
 399     fTrack[myIndex]->TOFsignal() = (Double_t)track->GetTOFsignal();
 400     for (k = 0; k < AliPID::kSPECIES; k++)
 401     {
 402       fTrack[myIndex]->TOFref(k)   = time[k];
 403       fTrack[myIndex]->TOFsigma(k) = (Double_t)fTOFmaker->GetExpectedSigma(fTrack[myIndex]->Prec(), time[k], AliPID::ParticleMass(k));
 404     }
 405   }
 406
 407   return kTRUE;
 408 }