PWGLF/SPECTRA/ParticleEfficiency/AliAnalysisTaskParticleEfficiency.cxx

   1 #include "AliAnalysisTaskParticleEfficiency.h"
   2 #include "AliESDEvent.h"
   3 #include "AliMCEvent.h"
   4 #include "AliStack.h"
   5 #include "AliInputEventHandler.h"
   6 #include "AliAnalysisManager.h"
   7 #include "AliESDVertex.h"
   8 #include "AliCentrality.h"
   9 #include "AliESDtrack.h"
  10 #include "AliESDtrackCuts.h"
  11 #include "TParticle.h"
  12 #include "TObjArray.h"
  13 #include "TList.h"
  14 #include "TH2F.h"
  15 #include "TDatabasePDG.h"
  16 #include "TParticlePDG.h"
  17
  18 ClassImp(AliAnalysisTaskParticleEfficiency)
  19
  20 //_______________________________________________________
  21
  22 AliAnalysisTaskParticleEfficiency::AliAnalysisTaskParticleEfficiency(const Char_t *partName) :
  23   AliAnalysisTaskSE(partName),
  24   fParticlePdgCode(0),
  25   fTrackCuts(NULL)
  26 {
  27
  28   /*
  29    * default constructor
  30    */
  31
  32   /* set particle PDG code */
  33   TParticlePDG *ppdg = TDatabasePDG::Instance()->GetParticle(partName);
  34   if (ppdg) fParticlePdgCode = ppdg->PdgCode();
  35
  36   /* init track cuts */
  37   fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
  38   fTrackCuts->SetEtaRange(-0.8, 0.8);
  39   fTrackCuts->SetPtRange(0.15, 1.e10);
  40
  41   /* create output */
  42   fHistoList = new TList();
  43   fHistoList->SetOwner(kTRUE);
  44   DefineOutput(1, TList::Class());
  45 }
  46
  47 //_______________________________________________________
  48
  49 AliAnalysisTaskParticleEfficiency::~AliAnalysisTaskParticleEfficiency()
  50 {
  51
  52   /*
  53    * default destructor
  54    */
  55
  56 }
  57
  58 //_______________________________________________________
  59
  60 void
  61 AliAnalysisTaskParticleEfficiency::UserCreateOutputObjects()
  62 {
  63
  64   /*
  65    * user create output objects
  66    */
  67
  68   fHistoEvents = new TH2F("hHistoEvents", ";centrality percentile;track multiplicity", 20, 0., 100., 200, 0., 2000.);
  69   fHistoList->Add(fHistoEvents);
  70
  71   fHistoGenerated = new TH2F("hHistoGenerated", ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
  72   fHistoList->Add(fHistoGenerated);
  73
  74   fHistoReconstructed = new TH2F("hHistoReconstructed", ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
  75   fHistoList->Add(fHistoReconstructed);
  76
  77   for (Int_t idau = 0; idau < 2; idau++) {
  78
  79     fHistoGeneratedDaughter[idau] = new TH2F(Form("hHistoGeneratedDaughter_%d", idau), ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
  80     fHistoList->Add(fHistoGeneratedDaughter[idau]);
  81
  82     fHistoReconstructedDaughter[idau] = new TH2F(Form("hHistoReconstructedDaughter_%d", idau), ";centrality percentile;p_{T} (GeV/c)", 20, 0., 100., 200, 0., 10.);
  83     fHistoList->Add(fHistoReconstructedDaughter[idau]);
  84
  85   }
  86
  87   PostData(1, fHistoList);
  88 }
  89
  90 //_______________________________________________________
  91
  92 void
  93 AliAnalysisTaskParticleEfficiency::UserExec(Option_t *)
  94 {
  95
  96   /*
  97    * user exec
  98    */
  99
 100   /*** INIT ***/
 101
 102   /* get ESD event */
 103   AliESDEvent *esdEvent = dynamic_cast<AliESDEvent *>(InputEvent());
 104   if (!esdEvent) return;
 105   /* get MC event */
 106   AliMCEvent *mcEvent = dynamic_cast<AliMCEvent *>(MCEvent());
 107   if (!mcEvent) return;
 108   /* get stack */
 109   AliStack *mcStack = mcEvent->Stack();
 110   if (!mcStack) return;
 111
 112   /*** EVENT SELECTION ***/
 113
 114   /* collision candidate */
 115   if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return;
 116   /* vertex selection */
 117   const AliESDVertex *vertex = esdEvent->GetPrimaryVertexTracks();
 118   if (vertex->GetNContributors() < 1) {
 119     vertex = esdEvent->GetPrimaryVertexSPD();
 120     if (vertex->GetNContributors() < 1) return;
 121   }
 122   /* z-vertex cut */
 123   if (TMath::Abs(vertex->GetZ()) > 10.) return;
 124   /* centrality selection */
 125   AliCentrality *centrality = esdEvent->GetCentrality();
 126   if (centrality->GetQuality() != 0) return;
 127   Double_t cent = centrality->GetCentralityPercentileUnchecked("V0M");
 128
 129   /* fill event histos */
 130   fHistoEvents->Fill(cent, fTrackCuts->CountAcceptedTracks(esdEvent));
 131
 132   /*** RECONSTRUCTED TRACKS ***/
 133
 134   TObjArray recoParticleArray;
 135
 136   /* loop over ESD tracks */
 137   for (Int_t itrk = 0; itrk < esdEvent->GetNumberOfTracks(); itrk++) {
 138
 139     /* get track */
 140     AliESDtrack *track = esdEvent->GetTrack(itrk);
 141     if (!track) continue;
 142     /* check accept track */
 143     if (!fTrackCuts->AcceptTrack(track)) continue;
 144
 145     /* get coresponding MC particle */
 146     TParticle *particle = mcStack->Particle(TMath::Abs(track->GetLabel()));
 147     if (!particle) continue;
 148
 149     /* add to reconstructed particle array */
 150     recoParticleArray.Add(particle);
 151   }
 152
 153   /*** MONTECARLO PARTICLES ***/
 154
 155   /* loop over MC stack */
 156   for (Int_t ipart = 0; ipart < mcStack->GetNtrack(); ipart++) {
 157     /* get particle */
 158     TParticle *particle = mcStack->Particle(ipart);
 159     if (!particle) continue;
 160     /* check PDG */
 161     if (particle->GetPdgCode() != fParticlePdgCode) continue;
 162     /* check rapidity */
 163     if (TMath::Abs(particle->Y()) > 0.5) continue;
 164
 165     /* switch PDG */
 166     switch (fParticlePdgCode) {
 167
 168     case 211: case -211: /* pions */
 169     case 321: case -321: /* kaons */
 170     case 2212: case -2212: /* protons */
 171
 172       /* check physical primary */
 173       if (!mcStack->IsPhysicalPrimary(ipart)) continue;
 174
 175       /* fill particle histos */
 176       fHistoGenerated->Fill(cent, particle->Pt());
 177       if (recoParticleArray.Contains(particle))
 178         fHistoReconstructed->Fill(cent, particle->Pt());
 179       break;
 180
 181     case 333: /* phi */
 182     case 313: /* K*0 */
 183
 184       /* check number of daughters */
 185       if (particle->GetNDaughters() != 2) continue;
 186       /* check daughter charge */
 187       TParticle *daughter0 =  mcStack->Particle(particle->GetDaughter(0));
 188       TParticle *daughter1 =  mcStack->Particle(particle->GetDaughter(1));
 189       if (TDatabasePDG::Instance()->GetParticle(daughter0->GetPdgCode())->Charge() == 0.) continue;
 190       if (TDatabasePDG::Instance()->GetParticle(daughter1->GetPdgCode())->Charge() == 0.) continue;
 191
 192       /* fill daughter histos */
 193       fHistoGeneratedDaughter[0]->Fill(cent, daughter0->Pt());
 194       if (recoParticleArray.Contains(daughter0))
 195         fHistoReconstructedDaughter[0]->Fill(cent, daughter0->Pt());
 196       fHistoGeneratedDaughter[1]->Fill(cent, daughter1->Pt());
 197       if (recoParticleArray.Contains(daughter1))
 198         fHistoReconstructedDaughter[1]->Fill(cent, daughter1->Pt());
 199
 200       /* fill particle histos */
 201       fHistoGenerated->Fill(cent, particle->Pt());
 202       if (recoParticleArray.Contains(daughter0) &&
 203           recoParticleArray.Contains(daughter1))
 204         fHistoReconstructed->Fill(cent, particle->Pt());
 205       break;
 206     }
 207
 208   }
 209
 210   PostData(1, fHistoList);
 211 }
 212