a24052f355cc143162160d702b39d1ea49b1c423
[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / extra / AliAnalysisTaskResonanceQA.cxx
1 #include <Riostream.h>
2
3 #include <TChain.h>
4 #include <TH1.h>
5 #include <TH2.h>
6 #include <TH3.h>
7 #include <TDatabasePDG.h>
8 #include <TParticlePDG.h>
9 #include <TParticle.h>
10 #include <TLorentzVector.h>
11
12 #include "AliESDpid.h"
13 #include "AliESDtrack.h"
14 #include "AliESDEvent.h"
15 #include "AliESDVertex.h"
16 #include "AliESDtrackCuts.h"
17
18 #include "AliPID.h"
19 #include "AliAnalysisTask.h"
20 #include "AliAnalysisManager.h"
21 #include "AliInputEventHandler.h"
22
23 #include "AliMCEvent.h"
24 #include "AliStack.h"
25 #include "AliMCParticle.h"
26 #include "AliGenEventHeader.h"
27
28 #include "AliAnalysisTaskResonanceQA.h"
29
30 //
31 // analysis task creating basic QA plots for resonance particles
32 // Author: Ayben Karasu Uysal
33 // Computes some QA histograms useful for checking productions and data
34 // both for counting resonances inside MC, and for checking PID performances
35 //
36 // Revised by A. Pulvirenti
37 //
38
39
40 ClassImp(AliAnalysisTaskResonanceQA)
41
42 //________________________________________________________________________
43 AliAnalysisTaskResonanceQA::AliAnalysisTaskResonanceQA(const char *name) :
44    AliAnalysisTaskSE(name), 
45    fT0(AliESDpid::kTOF_T0),
46    fPrimaryThr(1E-6),
47    fVz(10.0),
48    fOutputList(0),
49    fSelectedEvts(0),
50    fEventVz(0),
51    fdEdxTPC(0),
52    fdEdxITS(0),
53    fTOFpid(0),
54    fDCAXYvsPtBeforeCuts(0), 
55    fDCAZvsPtBeforeCuts(0),
56    fNClusterPtBeforeCuts(0), 
57    fNFindableClusterPtBeforeCuts(0), 
58    fNCrossedRowsTPCPtBeforeCuts(0),
59    fProducedParticles(0), 
60    fESD(0), 
61    fESDpid(0),      
62    fTrackCuts(0)
63 {
64 //   
65 // Constructor
66 //
67
68    // Define input and output slots here
69    // Input slot #0 works with a TChain
70    DefineInput(0, TChain::Class());
71    
72    // Output slot #0 id reserved by the base class for AOD
73    // Output slot #1 writes into a TH1 container
74    DefineOutput(1, TList::Class());
75    
76    // setup to NULL all remaining histos
77    Int_t i;
78    for (i = 0; i < kResonances; i++) fRsnYPt[0][i] = fRsnYPt[1][i] = 0;
79 }
80
81 //________________________________________________________________________
82 const char* AliAnalysisTaskResonanceQA::EvtName(EEvtType type) const
83 {
84 //
85 // Return a short string with name of resonance
86 //
87
88    switch(type) {
89       case kBadNoVertex   : return "NoVertex";
90       case kBadVertexOutZ : return "FarVertex";
91       case kGoodEvent     : return "Good";
92       default             : return "Unknown";
93    }
94 }
95
96 //________________________________________________________________________
97 const char* AliAnalysisTaskResonanceQA::RsnName(ERsn type) const
98 {
99 //
100 // Return a short string with name of resonance
101 //
102
103    switch(type) {
104       case kPhi       : return "Phi";
105       case kKStar0    : return "KStar0";
106       case kRho       : return "Rho";
107       case kLambdaStar: return "LambdaStar";
108       case kXiStar0   : return "XiStar0";
109       case kSigmaStarP: return "SigmaStarP";
110       case kSigmaStarM: return "SigmaStarM";
111       case kDeltaPP   : return "DeltaPP";
112       default         : return "Unknown";
113    }
114 }
115
116 //________________________________________________________________________
117 const char* AliAnalysisTaskResonanceQA::RsnSymbol(ERsn type) const
118 {
119 //
120 // Return a short string with name of resonance
121 //
122
123    switch(type) {
124       case kPhi       : return "#phi";
125       case kKStar0    : return "K^{*0}";
126       case kRho       : return "#rho";
127       case kLambdaStar: return "#Lambda^{*}";
128       case kXiStar0   : return "#Xi^{*0}";
129       case kSigmaStarP: return "#Sigma^{*+}";
130       case kSigmaStarM: return "#Sigma^{*-}";
131       case kDeltaPP   : return "#Delta^{++}";
132       default         : return "Unknown";
133    }
134 }
135
136 //________________________________________________________________________
137 Int_t AliAnalysisTaskResonanceQA::RsnPDG(ERsn type) const
138 {
139 //
140 // Return PDG code of resonance
141 //
142
143    switch(type) {
144       case kPhi       : return  333;
145       case kKStar0    : return  313;
146       case kRho       : return  113;
147       case kLambdaStar: return 3124;
148       case kXiStar0   : return 3324;
149       case kSigmaStarP: return 3224;
150       case kSigmaStarM: return 3114;
151       case kDeltaPP   : return 2224;
152       default         : return    0;
153    }
154 }
155
156 //________________________________________________________________________
157 void AliAnalysisTaskResonanceQA::UserCreateOutputObjects()
158 {
159 //
160 // Create histograms (called once)
161 //
162
163    Int_t i;
164    Int_t vESDdEdxBin = 1000;  Double_t vESDdEdxMin = 0;  Double_t vESDdEdxMax = 1000;
165    Int_t vESDQPBin = 900;     Double_t vESDQPMin = -4.5; Double_t vESDQPMax = 4.5;
166    Int_t vPtBin = 500;        Double_t vPtMin = 0.;      Double_t vPtMax = 10.;
167    Int_t vYBin = 600;         Double_t vYMin = -15.;     Double_t vYMax = 15.;
168
169    // standard objects
170    fESDpid = new AliESDpid();
171    fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
172    fTrackCuts->SetPtRange(0.1, 1E20);
173    fTrackCuts->SetEtaRange(-0.9, 0.9);
174
175    // create output list
176    fOutputList = new TList();
177    fOutputList->SetOwner();
178
179    // output histograms for PID signal
180    fSelectedEvts = new TH1I("fSelectedEvts", "fSelectedEvts;fSelectedEvts", 4, 0, 4);
181    fEventVz = new TH1F("fEventVz", "Vz distribution", 800, -40.0, 40.0);
182    fdEdxTPC = new TH2F("fdEdxTPC", "dEdxTPC, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
183    fdEdxITS = new TH2F("fdEdxITS", "dEdxITS, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
184    fTOFpid = new TH2F("fTOFpid", "TOF PID;p (GeV/c);#beta;", vESDQPBin, vESDQPMin, vESDQPMax, 300, 0.1, 1.2);
185    
186    // output histograms for track quality
187    fDCAXYvsPtBeforeCuts = new TH2F("DCAXYvsPtBeforeCuts", "DCAXYvsPtBeforeCuts;p_{t} (GeV/c);DCAXY", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
188    fDCAZvsPtBeforeCuts = new TH2F("DCAZvsPtBeforeCuts", "DCAZvsPtBeforeCuts;p_{t} (GeV/c);DCAZ", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
189    fNClusterPtBeforeCuts = new TH2F("NClusterPtBeforeCuts", "NClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
190    fNFindableClusterPtBeforeCuts = new TH2F("NFindableClusterPtBeforeCuts", "NFindableClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
191    fNCrossedRowsTPCPtBeforeCuts = new TH2F("NCrossedRowsTPCPtBeforeCuts", "NCrossedRowsTPCPtBeforeCuts;p_{t} (GeV/c);NCrossedRowsTPC", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
192    
193    // output histograms for resonances
194    fProducedParticles = new TH1I("ProducedParticles", "ProducedParticles;Produced Particles", kResonances, 0, kResonances);
195    for (i = 0; i < kResonances; i++) {
196       fRsnYPt[0][i] = new TH3F(Form("%s_all" , RsnName(i)), Form("%s -- ALL;p_{t} (GeV/c);Rapidity"    , RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
197       fRsnYPt[1][i] = new TH3F(Form("%s_prim", RsnName(i)), Form("%s -- PRIMARY;p_{t} (GeV/c);Rapidity", RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
198       fProducedParticles->GetXaxis()->SetBinLabel(i + 1, RsnSymbol(i));
199    }
200    
201    // add histogram to list
202    fOutputList->Add(fSelectedEvts);
203    fOutputList->Add(fEventVz);
204    fOutputList->Add(fdEdxTPC);
205    fOutputList->Add(fdEdxITS);
206    fOutputList->Add(fTOFpid);
207    fOutputList->Add(fDCAXYvsPtBeforeCuts);
208    fOutputList->Add(fDCAZvsPtBeforeCuts);
209    fOutputList->Add(fNClusterPtBeforeCuts);
210    fOutputList->Add(fNFindableClusterPtBeforeCuts);
211    fOutputList->Add(fNCrossedRowsTPCPtBeforeCuts);
212    fOutputList->Add(fProducedParticles);
213    for (i = 0; i < kResonances; i++) {
214       fOutputList->Add(fRsnYPt[0][i]);
215       fOutputList->Add(fRsnYPt[1][i]);
216    }
217    
218    PostData(1, fOutputList);
219 }
220
221 //________________________________________________________________________
222 void AliAnalysisTaskResonanceQA::UserExec(Option_t *)
223 {
224 //
225 // Execute computations
226 //
227
228    // first bin of this histogram counts processed events
229    // second bin counts events passing physics selection
230    // all events not passing physics selections are skipped
231    fSelectedEvts->Fill(0.1);
232    Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
233    if (!isSelected) return;
234    
235    // count CINT1B selected events
236    fSelectedEvts->Fill(1.1);
237    
238    // try to cast input event to ESD and loop on tracks
239    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
240    
241    // check if event is accepted w.r. to all additional selection criteria
242    Bool_t hasVertex = kFALSE, goodVertex = kFALSE;
243    Double_t vz = 1E20, vzsign = 1E20;
244    if (fESD) {
245       // check primary vertex:
246       // we accept only tracks/SPD with at least 1 contributor
247       const AliESDVertex *v0 = fESD->GetPrimaryVertexTracks();
248       if (v0) {
249          if (v0->GetNContributors() > 0) {
250             vzsign = v0->GetZv();
251             vz = TMath::Abs(vzsign);
252             hasVertex = kTRUE;
253          } else {
254             v0 = fESD->GetPrimaryVertexSPD();
255             if (v0) {
256                if (v0->GetNContributors() > 0) {
257                   vzsign = v0->GetZv();
258                   vz = TMath::Abs(vzsign);
259                   hasVertex = kTRUE;
260                }
261             }
262          }
263       }
264       if (hasVertex) {
265          goodVertex = (vz <= fVz);
266          fSelectedEvts->Fill(2.1);
267          fEventVz->Fill(vzsign);
268       }
269    }
270    
271    // determine event quality
272    Bool_t eventAccepted = kFALSE;
273    Int_t evtQuality = kEvtTypes;
274    if (hasVertex) {
275       if (goodVertex) {
276          evtQuality = kGoodEvent;
277          eventAccepted = kTRUE;
278       } else {
279          evtQuality = kBadVertexOutZ;
280       }
281    } else {
282       evtQuality = kBadNoVertex;
283    }
284    
285    // Message
286    static Int_t evNum = -1;
287    evNum++;
288    switch (evtQuality) {
289       case kBadNoVertex  : AliInfo(Form("Rejecting event #%5d because it has not a good vertex", evNum)); break;
290       case kBadVertexOutZ: AliInfo(Form("Rejecting event #%5d because it is out of range: vz = %5.3f vs. %5.3f", evNum, vz, fVz)); break;
291       case kGoodEvent    : AliDebugClass(1, Form("Accepting event #%5d", evNum)); break;
292       default            : AliError("This should never appear!!!");
293    }
294    
295    // if event is OK, loop on tracks to fill QA
296    if (fESD && eventAccepted) {
297       // if we arrive here, the event was processed successfully
298       // then we fill last bin in first histogram
299       fSelectedEvts->Fill(3.1);
300       
301       // settings for TOF time zero
302       if (fESD->GetTOFHeader()) 
303          fESDpid->SetTOFResponse(fESD, AliESDpid::kTOF_T0);
304       else {
305          Float_t t0spread[10];
306          Float_t intrinsicTOFres = 100;
307          for (Int_t i = 0; i < 10; i++) t0spread[i] = (TMath::Sqrt(fESD->GetSigma2DiamondZ())) / 0.03; //0.03 to convert from cm to ps
308          fESDpid->GetTOFResponse().SetT0resolution(t0spread);
309          fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
310          fESDpid->SetTOFResponse(fESD, fT0);
311       }
312       fESDpid->MakePID(fESD, kFALSE, 0.);
313
314       // use some variables to store useful info
315       /*MISC   */ Int_t    iTrack, nTracks = fESD->GetNumberOfTracks();
316       /*MISC   */ Double_t pt;
317       /*QUALITY*/ Float_t  impactXY, impactZ, nCrossedRowsTPC;
318       /*ITS    */ Double_t itsSignal = 0;
319       /*TPC    */ Double_t ptotSE = 0, tpcSignal, signSE;
320       /*TOF    */ Double_t momentum, length, time, beta; //, times[AliPID::kSPECIES], tzeroTrack = 0;
321
322       // loop on tracks
323       for (iTrack = 0; iTrack < nTracks; iTrack++) {
324          
325          // get track
326          AliESDtrack* track = fESD->GetTrack(iTrack);
327          if (!track) continue;
328          
329          // get transverse momentum (used for quality histograms)
330          pt = track->Pt();
331          
332          // get charge and reject neutrals
333          signSE = track->GetSign();
334          if (signSE == 0) continue;
335
336          // get DCA and TPC-related quality params
337          track->GetImpactParameters(impactXY, impactZ);
338          nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);
339
340          // fill quality histograms when status flags are OK for TPC in+refit and ITS refit
341          if (track->IsOn(AliESDtrack::kTPCin) && track->IsOn(AliESDtrack::kTPCrefit) && track->IsOn(AliESDtrack::kITSrefit)) {
342             fNClusterPtBeforeCuts        ->Fill(pt, track->GetTPCNcls());
343             fNFindableClusterPtBeforeCuts->Fill(pt, track->GetTPCNclsF());
344             fNCrossedRowsTPCPtBeforeCuts ->Fill(pt, nCrossedRowsTPC);
345             fDCAXYvsPtBeforeCuts         ->Fill(pt, impactXY);
346             fDCAZvsPtBeforeCuts          ->Fill(pt, impactZ);
347          }
348          
349          // from now on, work only with tracks passing standard cuts for 2010
350          if (!fTrackCuts->AcceptTrack(track)) continue;
351
352          // get momentum (used for PID histograms)
353          momentum = track->P();
354
355          // ITS
356          itsSignal = track->GetITSsignal();
357          fdEdxITS->Fill(signSE * momentum, itsSignal);
358          
359          // TPC
360          if (track->GetInnerParam()) {
361             AliExternalTrackParam trackIn1(*track->GetInnerParam());
362             ptotSE = trackIn1.GetP();
363             tpcSignal = track->GetTPCsignal();
364             fdEdxTPC->Fill(signSE * ptotSE, tpcSignal);
365          } // end TPC
366          
367          // TOF (requires checking some more status flags
368          if (track->IsOn(AliESDtrack::kTOFpid) && track->IsOn(AliESDtrack::kTOFout) && track->IsOn(AliESDtrack::kTIME) && !track->IsOn(AliESDtrack::kTOFmismatch)) {
369             length = track->GetIntegratedLength();
370             
371             // reject suspiciously small lengths/times
372             if (track->GetIntegratedLength() < 350.) continue;
373             if (track->GetTOFsignal() < 1E-6) continue;
374
375             // thhese are maybe not needed
376             // track->GetIntegratedTimes(times);
377             // tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
378             
379             // get TOF measured time and compute beta
380             time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(momentum);
381             beta = length / (2.99792457999999984e-02 * time);
382             fTOFpid->Fill(signSE * momentum, beta);
383          } // end TOF
384       } // end track loop
385    } // end ESD check
386    
387    // loop on MC, if MC event and stack are available
388    AliMCEvent *mcEvent = MCEvent();
389    if (mcEvent) {
390       // MC primary vertex
391       TArrayF mcVertex(3);
392       AliGenEventHeader *genEH = mcEvent->GenEventHeader();
393       genEH->PrimaryVertex(mcVertex);
394       // loop on particles
395       AliStack *stack = mcEvent->Stack();
396       if (stack) {
397          Int_t iTrack, irsn, pdg, mother, motherPDG, nTracks = stack->GetNtrack();
398          Double_t dx, dy, dz, dist;
399          TLorentzVector vprod;
400          for (iTrack = 0; iTrack < nTracks; iTrack++) {
401             
402             // get particle
403             TParticle *part = stack->Particle(iTrack);
404             if (!part) {
405                AliError(Form("Could not receive track %d", iTrack));
406                continue;
407             }
408             
409             // get PDG code of particle and check physical primary
410             pdg = part->GetPdgCode();
411             mother = part->GetFirstMother();
412             motherPDG = 0;
413                         
414             // loop over possible resonances
415             for (irsn = 0; irsn < kResonances; irsn++) {
416                if (TMath::Abs(pdg) == RsnPDG(irsn)) {
417                   // debug check
418                   if (mother >= 0) {
419                      TParticle *p = stack->Particle(mother);
420                      motherPDG = p->GetPdgCode();
421                   }
422                   part->ProductionVertex(vprod);
423                   dx = mcVertex[0] - vprod.X();
424                   dy = mcVertex[1] - vprod.Y();
425                   dz = mcVertex[2] - vprod.Z();
426                   dist = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
427                   AliDebugClass(1, Form("PDG = %d -- mother ID = %d, pdg = %d -- dist. from MC vertex = %f -- prod. time = %f", pdg, mother, motherPDG, dist, vprod.T())); 
428                   // global counter is always filled
429                   fRsnYPt[0][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
430                   // counter for primaries is filled only if mother is less than 0 (primary)
431                   if (dist < fPrimaryThr) {
432                      fRsnYPt[1][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
433                   }
434                   // fill the global histogram
435                   fProducedParticles->Fill(irsn + 1);
436                   // if one PDG match was found, no need to continue loop
437                   break;
438                }
439             }
440          }
441       }
442    }
443    
444    PostData(1, fOutputList);
445 }
446
447 //________________________________________________________________________
448 void AliAnalysisTaskResonanceQA::Terminate(Option_t *)
449 {
450 //
451 // Termination
452 // Quickly check that the output list is there
453 //
454   
455    fOutputList = dynamic_cast<TList*>(GetOutputData(1));
456    if (!fOutputList) {
457       AliError("Output list not found!!!");
458       return;
459    }
460 }