PWGLF/SPECTRA/PiKaPr/TOF/PbPb276/task/AliAnalysisTaskTOFSpectraPbPb.cxx

   1 #include "AliAnalysisTaskTOFSpectraPbPb.h"
   2 #include "AliESDEvent.h"
   3 #include "AliMCEvent.h"
   4 #include "AliStack.h"
   5 #include "AliPhysicsSelection.h"
   6 #include "AliESDtrackCuts.h"
   7 #include "AliESDpid.h"
   8 #include "AliTOFcalib.h"
   9 #include "AliTOFT0maker.h"
  10 #include "TList.h"
  11 #include "TH1F.h"
  12 #include "TH2F.h"
  13 #include "AliCDBManager.h"
  14 #include "AliLog.h"
  15 #include "AliESDtrack.h"
  16 #include "TObjArray.h"
  17 #include "TLorentzVector.h"
  18 #include "TParticle.h"
  19 #include "TDatabasePDG.h"
  20 #include "TParticlePDG.h"
  21 #include "AliTOFRunParams.h"
  22 #include "AliAnalysisTrack.h"
  23 #include "AliAnalysisParticle.h"
  24 #include "AliAnalysisEvent.h"
  25 #include "TClonesArray.h"
  26 #include "AliAnalysisManager.h"
  27 #include "AliAODHandler.h"
  28 #include "TTree.h"
  29 #include "AliCentrality.h"
  30 #include "AliInputEventHandler.h"
  31 #include "AliVEvent.h"
  32
  33 ClassImp(AliAnalysisTaskTOFSpectraPbPb)
  34
  35 //_______________________________________________________
  36
  37 AliAnalysisTaskTOFSpectraPbPb::AliAnalysisTaskTOFSpectraPbPb() :
  38   AliAnalysisTaskSE("TOFSpectraPbPb"),
  39   fInitFlag(kFALSE),
  40   fMCFlag(kFALSE),
  41   fMCTuneFlag(kFALSE),
  42   fPbPbFlag(kFALSE),
  43   fVertexSelectionFlag(kFALSE),
  44   fRunNumber(0),
  45   fESDEvent(NULL),
  46   fMCEvent(NULL),
  47   fMCStack(NULL),
  48   fTrackCuts(NULL),
  49   fESDpid(new AliESDpid()),
  50   fIsCollisionCandidate(kFALSE),
  51   fHasVertex(kFALSE),
  52   fVertexZ(0.),
  53   fMCTimeZero(0.),
  54   fCentrality(NULL),
  55   fAnalysisEvent(new AliAnalysisEvent()),
  56   fAnalysisTrackArray(new TClonesArray("AliAnalysisTrack")),
  57   fAnalysisTrack(new AliAnalysisTrack()),
  58   fAnalysisParticleArray(new TClonesArray("AliAnalysisParticle")),
  59   fAnalysisParticle(new AliAnalysisParticle()),
  60   fTOFcalib(new AliTOFcalib()),
  61   fTOFT0maker(new AliTOFT0maker(fESDpid)),
  62   fTimeResolution(80.),
  63   fVertexCut(10.0),
  64   fRapidityCut(0.5),
  65   fHistoList(new TList()),
  66   fMCHistoList(new TList())
  67 {
  68   /*
  69    * default constructor
  70    */
  71
  72 }
  73
  74 //_______________________________________________________
  75
  76 AliAnalysisTaskTOFSpectraPbPb::~AliAnalysisTaskTOFSpectraPbPb()
  77 {
  78   /*
  79    * default destructor
  80    */
  81
  82   if (fTrackCuts) delete fTrackCuts;
  83   delete fESDpid;
  84   delete fTOFcalib;
  85   delete fTOFT0maker;
  86   delete fHistoList;
  87   delete fMCHistoList;
  88 }
  89
  90 //_______________________________________________________
  91
  92 void
  93 AliAnalysisTaskTOFSpectraPbPb::UserCreateOutputObjects()
  94 {
  95   /*
  96    * user create output objects
  97    */
  98
  99   /* output tree */
 100   OutputTree()->Branch("AnalysisEvent", "AliAnalysisEvent", &fAnalysisEvent);
 101   OutputTree()->Branch("AnalysisTrack", "TClonesArray", &fAnalysisTrackArray);
 102   if (fMCFlag)
 103     OutputTree()->Branch("AnalysisParticle", "TClonesArray", &fAnalysisParticleArray);
 104 }
 105
 106 //_______________________________________________________
 107
 108 Bool_t
 109 AliAnalysisTaskTOFSpectraPbPb::InitRun()
 110 {
 111   /*
 112    * init run
 113    */
 114
 115   /* get ESD event */
 116   fESDEvent = dynamic_cast<AliESDEvent *>(InputEvent());
 117   if (!fESDEvent) {
 118     AliError("cannot get ESD event");
 119     return kFALSE;
 120   }
 121   /* get run number */
 122   Int_t runNb = fESDEvent->GetRunNumber();
 123   /* check run already initialized */
 124   if (fInitFlag && fRunNumber == runNb) return kTRUE;
 125   /* init cdb */
 126   AliCDBManager *cdb = AliCDBManager::Instance();
 127   cdb->SetDefaultStorage("raw://");
 128   cdb->SetRun(runNb);
 129   /* init TOF calib */
 130   if (!fTOFcalib->Init()) {
 131     AliError("cannot init TOF calib");
 132     return kFALSE;
 133   }
 134   AliInfo(Form("initialized for run %d", runNb));
 135   fInitFlag = kTRUE;
 136   fRunNumber = runNb;
 137   return kTRUE;
 138 }
 139
 140 //_______________________________________________________
 141
 142 Bool_t
 143 AliAnalysisTaskTOFSpectraPbPb::InitEvent()
 144 {
 145   /*
 146    * init event
 147    */
 148
 149   /* get ESD event */
 150   fESDEvent = dynamic_cast<AliESDEvent *>(InputEvent());
 151   if (!fESDEvent) return kFALSE;
 152   /* get MC event */
 153   if (fMCFlag) {
 154     fMCEvent = dynamic_cast<AliMCEvent *>(MCEvent());
 155     if (!fMCEvent) return kFALSE;
 156   }
 157   /* get stack */
 158   if (fMCFlag) {
 159     fMCStack = fMCEvent->Stack();
 160     if (!fMCStack) return kFALSE;
 161   }
 162   /* event selection */
 163   fIsCollisionCandidate = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 164   /* vertex selection */
 165   const AliESDVertex *vertex = fESDEvent->GetPrimaryVertexTracks();
 166   if (vertex->GetNContributors() < 1) {
 167     vertex = fESDEvent->GetPrimaryVertexSPD();
 168     if (vertex->GetNContributors() < 1) fHasVertex = kFALSE;
 169     else fHasVertex = kTRUE;
 170   }
 171   else fHasVertex = kTRUE;
 172   if (fVertexSelectionFlag && (!fHasVertex || TMath::Abs(vertex->GetZ()) > fVertexCut)) return kFALSE;
 173   fVertexZ = vertex->GetZ();
 174   /* centrality in PbPb */
 175   if (fPbPbFlag) {
 176     fCentrality = fESDEvent->GetCentrality();
 177   }
 178
 179   /* calibrate ESD (also in MC to correct TExp) */
 180   fTOFcalib->CalibrateESD(fESDEvent);
 181
 182   /* check MC flags */
 183   if (fMCFlag && fMCTuneFlag)
 184     fMCTimeZero = fTOFcalib->TuneForMC(fESDEvent, fTimeResolution);
 185
 186 #if 0 /* NEW METHOD */
 187   /* compute TOF-T0, fill ESD and make PID */
 188   fTOFT0maker->ComputeT0TOF(fESDEvent);
 189   fTOFT0maker->WriteInESD(fESDEvent);
 190   fESDpid->SetTOFResponse(fESDEvent, AliESDpid::kTOF_T0);
 191   fESDpid->MakePID(fESDEvent, kFALSE, 0.);
 192 #endif
 193
 194 #if 1 /* OLD METHOD */
 195   /* compute and apply TOF-T0 */
 196   fTOFT0maker->ComputeT0TOF(fESDEvent);
 197   fESDpid->MakePID(fESDEvent, kFALSE, 0.);
 198 #endif
 199
 200   return kTRUE;
 201 }
 202
 203 //_______________________________________________________
 204
 205 Bool_t
 206 AliAnalysisTaskTOFSpectraPbPb::HasPrimaryDCA(AliESDtrack *track)
 207 {
 208   /*
 209    * has primary DCA
 210    */
 211
 212   // cut on transverse impact parameter
 213   Float_t d0z0[2],covd0z0[3];
 214   track->GetImpactParameters(d0z0, covd0z0);
 215   Float_t sigma= 0.0050 + 0.0060 / TMath::Power(track->Pt(), 0.9);
 216   Float_t d0max = 7. * sigma;
 217   //
 218   Float_t sigmaZ = 0.0146 + 0.0070 / TMath::Power(track->Pt(), 1.114758);
 219   if (track->Pt() > 1.) sigmaZ = 0.0216;
 220   Float_t d0maxZ = 5. * sigmaZ;
 221   //
 222   if(TMath::Abs(d0z0[0]) > d0max || TMath::Abs(d0z0[1]) > d0maxZ)
 223     return kFALSE;
 224
 225   /* primary DCA ok */
 226   return kTRUE;
 227 }
 228
 229 //_______________________________________________________
 230
 231 void
 232 AliAnalysisTaskTOFSpectraPbPb::UserExec(Option_t *option)
 233 {
 234   /*
 235    * user exec
 236    */
 237
 238   /*** INITIALIZATION ***/
 239
 240   /* unset fill AOD */
 241   ((AliAODHandler*)(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()))->SetFillAOD(kFALSE);
 242
 243   /* init run */
 244   if (!InitRun()) return;
 245   /* init event */
 246   if (!InitEvent()) return;
 247
 248   /* set fill AOD */
 249   ((AliAODHandler*)(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()))->SetFillAOD(kTRUE);
 250
 251   /*** MC PRIMARY PARTICLES ***/
 252
 253   if (fMCFlag) {
 254
 255     /* reset track array */
 256     fAnalysisParticleArray->Clear();
 257
 258     /* loop over primary particles */
 259     Int_t nPrimaries = fMCStack->GetNprimary();
 260     TParticle *particle;
 261     TParticlePDG *particlePDG;
 262     /* loop over primary particles */
 263     for (Int_t ipart = 0; ipart < nPrimaries; ipart++) {
 264       /* check primary */
 265       if (!fMCStack->IsPhysicalPrimary(ipart)) continue;
 266       /* get particle */
 267       particle = fMCStack->Particle(ipart);
 268       if (!particle) continue;
 269       /* check rapidity and pt cuts */
 270       if (TMath::Abs(particle->Y()) > 0.5) continue;
 271       if (particle->Pt() < 0.15) continue;
 272       /* get particlePDG */
 273       particlePDG = particle->GetPDG();
 274       if (!particlePDG) continue;
 275       /* check charged */
 276       if (particlePDG->Charge() == 0.) continue;
 277       /* update and add analysis particle */
 278       fAnalysisParticle->Update(particle, ipart);
 279       new ((*fAnalysisParticleArray)[fAnalysisParticleArray->GetEntries()]) AliAnalysisParticle(*fAnalysisParticle);
 280     } /* end of loop over primary particles */
 281   }
 282
 283   /*** GLOBAL EVENT INFORMATION ***/
 284
 285   fAnalysisEvent->Reset();
 286   /* update global event info */
 287   fAnalysisEvent->SetIsCollisionCandidate(fIsCollisionCandidate);
 288   fAnalysisEvent->SetHasVertex(fHasVertex);
 289   fAnalysisEvent->SetVertexZ(fVertexZ);
 290   fAnalysisEvent->SetMCTimeZero(fMCTimeZero);
 291   /* update TOF event info */
 292   for (Int_t i = 0; i < 10; i++) {
 293     fAnalysisEvent->SetTimeZeroTOF(i, fESDpid->GetTOFResponse().GetT0bin(i));
 294     fAnalysisEvent->SetTimeZeroTOFSigma(i, fESDpid->GetTOFResponse().GetT0binRes(i));
 295   }
 296   /* update T0 event info */
 297   for (Int_t i = 0; i < 3; i++)
 298     fAnalysisEvent->SetTimeZeroT0(i, fESDEvent->GetT0TOF(i));
 299
 300   /* update centrality info in PbPb */
 301   if (fPbPbFlag) {
 302     fAnalysisEvent->SetCentralityQuality(fCentrality->GetQuality());
 303     for (Int_t icent = 0; icent < AliAnalysisEvent::kNCentralityEstimators; icent++)
 304       fAnalysisEvent->SetCentralityPercentile(icent, fCentrality->GetCentralityPercentileUnchecked(AliAnalysisEvent::fgkCentralityEstimatorName[icent]));
 305   }
 306
 307   /*** RECONSTRUCTED TRACKS ***/
 308
 309   /* reset track array */
 310   fAnalysisTrackArray->Clear();
 311
 312   /* loop over ESD tracks */
 313   Int_t nTracks = fESDEvent->GetNumberOfTracks();
 314   AliESDtrack *track;
 315   for (Int_t itrk = 0; itrk < nTracks; itrk++) {
 316     /* get track */
 317     track = fESDEvent->GetTrack(itrk);
 318     if (!track) continue;
 319     /* check accept track */
 320     if (!fTrackCuts->AcceptTrack(track)) continue;
 321     /* update and add analysis track */
 322     fAnalysisTrack->Update(track, fMCStack, fMCEvent);
 323     new ((*fAnalysisTrackArray)[fAnalysisTrackArray->GetEntries()]) AliAnalysisTrack(*fAnalysisTrack);
 324
 325   } /* end of loop over ESD tracks */
 326
 327 }
 328