Cleanup the code. Fix memory leak. Now inherit from AliAnalysisTaskSE (Antoine, Phili...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliAnalysisTaskCheckV0tender.cxx
1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  *
14 **************************************************************************/
15
16 /* $Id$ */
17
18 //
19 // Task fir checking the performance of the V0 tender
20 // 
21 // 
22 // Authors
23 //   Matus Kalisky <matus.kalisky@cern.ch>
24 //
25
26 #include <sstream>
27
28 #include <TH1F.h>
29 #include <TList.h>
30
31 #include "AliAnalysisManager.h"
32 #include "AliMCEventHandler.h"
33 #include "AliESDInputHandler.h"
34 #include "AliESDv0.h"
35 #include "AliESDtrack.h"
36 #include "AliMCParticle.h"
37 #include "AliMCEvent.h"
38
39 #include "AliHFEtools.h"
40 #include "AliHFEcollection.h"
41
42 #include "AliAnalysisTaskCheckV0tender.h"
43
44 ClassImp(AliAnalysisTaskCheckV0tender)
45
46 //__________________________________________________________
47 AliAnalysisTaskCheckV0tender::AliAnalysisTaskCheckV0tender():
48   AliAnalysisTaskSE("CheckV0tenderTask")
49   , fOutput(0x0)
50   , fColl(0x0)
51   , fCollMC(0x0)
52   , fEvents(0x0)
53 {
54   //
55   // Default Constructor
56   //
57 }
58 //__________________________________________________________
59 AliAnalysisTaskCheckV0tender::AliAnalysisTaskCheckV0tender(const Char_t *name):
60   AliAnalysisTaskSE(name)
61   , fOutput(0x0)
62   , fColl(0x0)
63   , fCollMC(0x0)
64   , fEvents(0x0)
65 {
66   //
67   // Default Constructor
68   //
69   DefineOutput(1, TH1F::Class());
70   DefineOutput(2, TList::Class());
71
72 }
73 //__________________________________________________________
74 AliAnalysisTaskCheckV0tender::~AliAnalysisTaskCheckV0tender(){
75   //
76   // Destructor
77   //
78
79   if (fOutput) delete fOutput;
80   //if (fColl) delete fColl;
81   //if (fCollMC) delete fCollMC;
82   if (fEvents) delete fEvents;
83 }
84 //__________________________________________________________
85 void AliAnalysisTaskCheckV0tender::UserCreateOutputObjects(){
86   //
87   // prepare output objects
88   //
89
90   fOutput = new TList();
91   fOutput->SetOwner();
92   // Counter for number of events
93   fEvents = new TH1I("nEvents", "NumberOfEvents", 1, 1, 2);
94
95   fColl = new AliHFEcollection("V0_QA", "tender V0s for data");
96   fCollMC = new AliHFEcollection("V0_MC_QA", "tender V0s for MC");
97
98   // 
99   // Data histos
100   //
101   
102   // total number of tagged V0s
103   fColl->CreateTH1F("h_NumberOf_V0s", "Number of tagged V0s; type; counts", 4, -0.5, 3.5);
104   // pT spectra of the tagged V0s
105   fColl->CreateTH1F("h_Gamma_pt", "p_{T} spectrum of tagged gammas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
106   fColl->CreateTH1F("h_K0_pt", "p_{T} spectrum of tagged K0s; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
107   fColl->CreateTH1F("h_Lambda_pt", "p_{T} spectrum of tagged Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
108   fColl->CreateTH1F("h_ALambda_pt", "p_{T} spectrum of tagged A-Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
109   // invariant mass of the V0s
110   fColl->CreateTH1F("h_Gamma_mass", "Inv. Mass of the ", 100, 0., 0.1);
111   fColl->CreateTH1F("h_K0_mass", "Inv. Mass of the ", 100, 0.45, 0.55);
112   fColl->CreateTH1F("h_Lambda_mass", "Inv. Mass of the ", 100, 1.08, 1.14);
113   fColl->CreateTH1F("h_ALambda_mass", "Inv. Mass of the ", 100, 1.08, 1.14);
114
115   // total number of tagged daughter particles (should correlate with number of V0s !
116   fColl->CreateTH1F("h_NumberOfDaughters", "Number of tagged daughters; type; counts", 3, -0.5, 2.5);
117   // pT spectra of tagged daughter particles
118   fColl->CreateTH1F("h_Electron_pt", "p_{T} spectrum of tagged; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
119   fColl->CreateTH1F("h_Pion_pt", "p_{T} spectrum of tagged; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
120   fColl->CreateTH1F("h_Proton_pt", "p_{T} spectrum of tagged; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
121
122   //
123   // MC histos
124   //
125
126   // pT spectra of the tagged V0s
127   fCollMC->CreateTH1F("h_Gamma_pt_S", "MC-S: p_{T} spectrum of tagged gammas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
128   fCollMC->CreateTH1F("h_K0_pt_S", "MC-S: p_{T} spectrum of tagged K0s; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
129   fCollMC->CreateTH1F("h_Lambda_pt_S", "MC-S: p_{T} spectrum of tagged Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
130   fCollMC->CreateTH1F("h_ALambda_pt_S", "MC-S: p_{T} spectrum of tagged A-Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0); 
131   fCollMC->CreateTH1F("h_Gamma_pt_B", "MC-B: p_{T} spectrum of tagged gammas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
132   fCollMC->CreateTH1F("h_K0_pt_B", "MC-B: p_{T} spectrum of tagged K0s; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
133   fCollMC->CreateTH1F("h_Lambda_pt_B", "MC-B: p_{T} spectrum of tagged Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
134   fCollMC->CreateTH1F("h_ALambda_pt_B", "MC-B: p_{T} spectrum of tagged A-Lambdas; p_{T} (GeV/c); counts", 20, 0.1, 20, 0); 
135   // invariant mass of the V0s
136   fCollMC->CreateTH1F("h_Gamma_mass_S", "MC-S: Inv. Mass of the gamma; m (GeV/c^{2}); counts", 100, 0., 0.1);
137   fCollMC->CreateTH1F("h_K0_mass_S", "MC-S: Inv. Mass of the K0; m (GeV/c^{2}); counts", 100, 0.45, 0.55);
138   fCollMC->CreateTH1F("h_Lambda_mass_S", "MC-S: Inv. Mass of the Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
139   fCollMC->CreateTH1F("h_ALambda_mass_S", "MC-S: Inv. Mass of the A-Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
140   fCollMC->CreateTH1F("h_Gamma_mass_B", "MC-B: Inv. Mass of the gamma; m (GeV/c^{2}); counts", 100, 0., 0.1);
141   fCollMC->CreateTH1F("h_K0_mass_B", "MC-B: Inv. Mass of the K0; m (GeV/c^{2}); counts", 100, 0.45, 0.55);
142   fCollMC->CreateTH1F("h_Lambda_mass_B", "MC-B: Inv. Mass of the Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
143   fCollMC->CreateTH1F("h_ALambda_mass_B", "MC-B: Inv. Mass of the A-Lambda; m (GeV/c^{2}); counts", 100, 1.08, 1.14);
144   // pT spectra of tagged daughter particles
145   fCollMC->CreateTH1F("h_Electron_pt_S", "MC-S: p_{T} spectrum of tagged electrons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
146   fCollMC->CreateTH1F("h_Pion_pt_S", "MC-S: p_{T} spectrum of tagged pions; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
147   fCollMC->CreateTH1F("h_Proton_pt_S", "MC-S: p_{T} spectrum of tagged protons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
148   fCollMC->CreateTH1F("h_Electron_pt_B", "MC-B: p_{T} spectrum of tagged electrons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
149   fCollMC->CreateTH1F("h_Pion_pt_B", "MC-B: p_{T} spectrum of tagged pions; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
150   fCollMC->CreateTH1F("h_Proton_pt_B", "MC-B: p_{T} spectrum of tagged protons; p_{T} (GeV/c); counts", 20, 0.1, 20, 0);
151
152
153   TList *tmp = fColl->GetList();
154   tmp->SetName(fColl->GetName());
155   fOutput->Add(tmp);
156   tmp = 0x0;
157   tmp = fCollMC->GetList();
158   tmp->SetName(fCollMC->GetName());
159   fOutput->Add(tmp);
160   
161   
162 }
163 //__________________________________________________________
164 void AliAnalysisTaskCheckV0tender::UserExec(Option_t *){
165   //
166   // Event Loop
167   // 
168
169   AliMCEventHandler* mcHandler = (dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()));
170   //AliESDInputHandler *inh = dynamic_cast<AliESDInputHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
171   //AliESDpid *workingPID = NULL;
172   //if(inh && (workingPID = inh->GetESDpid())) workingPID = inh->GetESDpid();
173   //else workingPID = AliHFEtools::GetDefaultPID(mcHandler ? kTRUE : kFALSE);
174      
175   // check the MC data
176   if(fMCEvent && !mcHandler ) return;
177   if(fMCEvent &&  !mcHandler->InitOk() ) return;
178   if(fMCEvent &&  !mcHandler->TreeK() ) return;
179   if(fMCEvent &&  !mcHandler->TreeTR() ) return;
180
181   ProcessV0s();
182   ProcessDaughters();
183
184   if(fMCEvent){
185     ProcessV0sMC();
186     ProcessDaughtersMC();
187   }
188
189   fEvents->Fill(1.1);
190   PostData(1, fEvents);
191   PostData(2, fOutput);
192 }
193 //__________________________________________________________
194 void AliAnalysisTaskCheckV0tender::Terminate(Option_t *){
195   //
196   // Do Post Processing
197   //
198 }
199 //__________________________________________________________
200 void AliAnalysisTaskCheckV0tender::ProcessV0s(){
201   //
202   // loop over the V0s and extract the information about
203   // the V0s tagged by the V0 tender
204   //
205
206   const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
207   TString name;
208
209   Int_t nV0s = fInputEvent->GetNumberOfV0s();
210   for(Int_t i=0; i<nV0s; ++i){
211     AliESDv0 *esdV0 = (static_cast<AliESDEvent *>(fInputEvent))->GetV0(i);
212     if(!esdV0) continue;
213     if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
214     Int_t pid = GetTenderPidV0(esdV0);
215     if(pid < 0) continue;
216     fColl->Fill("h_NumberOf_V0s", pid);
217     Float_t pT = esdV0->Pt();
218     name = "h_" + type[pid] + "_pt";
219     fColl->Fill(name, pT);
220     Float_t mass = MassV0(esdV0, pid);    
221     name = "h_" + type[pid] + "_mass";
222     //printf(" -D: name: %s \n", name.Data());
223     fColl->Fill(name, mass);
224   }
225
226 }
227 //__________________________________________________________
228 void AliAnalysisTaskCheckV0tender::ProcessDaughters(){
229   //
230   // produce some check plots for V0 tender tagged single tracks
231   //
232
233   const TString type[3] = {"Electron", "Pion", "Proton"};
234   TString name;
235
236   Int_t nTracks = fInputEvent->GetNumberOfTracks();
237   for(Int_t i=0; i<nTracks; ++i){
238     AliESDtrack *track = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(i));
239     if(!track) continue;
240     Int_t  pid = GetTenderPidDaughter(track);
241     if(pid < 0) continue;
242     fColl->Fill("h_NumberOfDaughters", pid*1.0);
243     Float_t pT = track->Pt();
244     name = "h_" + type[pid] + "_pt";
245     fColl->Fill(name, pT);
246     
247   }
248 }
249 //__________________________________________________________
250 void AliAnalysisTaskCheckV0tender::ProcessV0sMC(){
251   //
252   // check all V0tender selected V0 on their true identity
253   // 
254
255   const Int_t pid2pdg[4] = {22, 310, 3122, -3122};
256   const TString type[4] = {"Gamma", "K0", "Lambda", "ALambda"};
257   TString name;
258   Int_t nV0s = fInputEvent->GetNumberOfV0s();
259
260   Int_t nTracks = fInputEvent->GetNumberOfTracks();
261
262   // V0 loop
263   for(Int_t i=0; i<nV0s; ++i){
264     Bool_t id = kFALSE;
265     AliESDv0 *esdV0 = (static_cast<AliESDEvent *>(fInputEvent))->GetV0(i);
266     if(!esdV0) continue;
267     if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
268     Int_t pid = GetTenderPidV0(esdV0);
269     if(pid < 0) continue;
270
271     // both ESD daughtr tracks
272     Int_t iN, iP;
273     iN = iP = -1;
274     iP = esdV0->GetPindex();
275     iN = esdV0->GetNindex();    
276     if(iN < 0 || iP < 0) continue;
277     if(iN >= nTracks || iP >= nTracks) continue;    
278     AliESDtrack *dP, *dN;
279     dP = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(iP));
280     dN = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(iN));  
281     if(!dN || !dP) continue;
282
283     // MC labels of the daughter tracks
284     Int_t lN, lP;
285     lN = dN->GetLabel();
286     lP = dP->GetLabel();
287     if(lN < 0 || lP < 0) continue;
288
289     // MC daughter particles
290     AliMCParticle *mcP, *mcN;
291     mcP = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lP));
292     mcN = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lN));
293     if(!mcP || !mcN) continue;
294
295     // labels of the mother particles
296     Int_t lPm, lNm;
297     lPm = mcP->GetMother();
298     lNm = mcN->GetMother();
299     if(lPm < 0) continue;
300     AliMCParticle *m = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lPm)); 
301     if(!m) continue;
302     Int_t pdg = m->PdgCode();
303     if((lPm == lNm) && (pdg == pid2pdg[pid])) id = kTRUE;
304
305     if(id) name = "h_" + type[pid] + "_pt_S"; 
306     else name = "h_" + type[pid] + "_pt_B"; 
307     Float_t pT = esdV0->Pt();
308     fCollMC->Fill(name, pT);
309
310     if(id) name = "h_" + type[pid] + "_mass_S";
311     else name = "h_" + type[pid] + "_mass_B";
312     Float_t mass = MassV0(esdV0, pid);
313     fCollMC->Fill(name, mass);
314
315    }
316 }
317 //__________________________________________________________
318 void AliAnalysisTaskCheckV0tender::ProcessDaughtersMC(){
319   //
320   // check the identity of the V0tender selected V0 daughters
321   // !!! for positive check only the true identity plays a role here, 
322   // not the true V0 mother identity (e.g. selected electron could come
323   // from primary vertex or pion dalitz deca or true gamma conversion) !!!
324   //
325
326   const Int_t pid2pdg [3] = {11, 211, 2212};
327   const TString type[3] = {"Electron", "Pion", "Proton"};
328   TString name;
329
330
331   Int_t nTracks = fInputEvent->GetNumberOfTracks();
332   for(Int_t i=0; i<nTracks; ++i){
333     Bool_t id = kFALSE;
334     AliESDtrack *track = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(i));
335     if(!track) continue;
336     Int_t  pid = GetTenderPidDaughter(track);
337     if(pid < 0) continue;
338     Float_t pT = track->Pt();
339     Int_t label = track->GetLabel();
340     if(label < 0) continue;
341     AliMCParticle *mcp = 0x0;
342     mcp = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(label));
343     if(!mcp) continue;
344     Int_t pdg = TMath::Abs(mcp->PdgCode());
345     if(pdg == pid2pdg[pid]) id = kTRUE;
346     if(id) name = "h_" + type[pid] + "_pt_S";
347     else name = "h_" + type[pid] + "_pt_B";
348     fCollMC->Fill(name, pT);
349   }
350 }
351 //__________________________________________________________
352 Int_t AliAnalysisTaskCheckV0tender::GetTenderPidV0(AliESDv0 * const v0){
353   //
354   // retrieve the PID nformation stored in the status flags by the train tender
355   // 
356   Int_t pid = -1;
357   if(!v0){
358     return pid;
359   }
360   Int_t nTimes = 0;
361   if(v0->TestBit(BIT(14))){
362     pid = 0;
363     nTimes++;
364   }
365   if(v0->TestBit(BIT(15))){
366     pid = 1;
367     nTimes++;
368   }
369   if(v0->TestBit(BIT(16))){
370     pid = 2;
371     nTimes++;
372   }
373   if(v0->TestBit(BIT(17))){
374     pid = 3;
375     nTimes++;
376   }
377   if(nTimes > 1){
378     AliWarning("V0 track labeled multiple times by the V0 tender");
379     pid = -1;
380   }    
381
382   //printf(" -D: pid: %i \n", pid);
383
384   return pid;
385 }
386 //__________________________________________________________
387 Int_t AliAnalysisTaskCheckV0tender::GetTenderPidDaughter(AliESDtrack * const track){
388   //
389   // retrieve the PID nformation stored in the status flags by the train tender
390   // 
391
392   Int_t pid = -1;
393   if(!track){
394     return pid;
395   }
396   Int_t nTimes = 0;
397   if(track->TestBit(BIT(14))){
398     pid = 0;
399     nTimes++;
400   }
401   if(track->TestBit(BIT(15))){
402     pid = 1;
403     nTimes++;
404   }
405   if(track->TestBit(BIT(16))){
406     pid = 2;
407     nTimes++;
408   }
409   if(nTimes > 1){
410     AliWarning("V0 track labeled multiple times by the V0 tender");
411     pid = -1;
412   }    
413   return pid;
414 }
415 //__________________________________________________________
416 Float_t AliAnalysisTaskCheckV0tender::MassV0(AliESDv0 * const v0, Int_t id){
417   //
418   // Get the V0 effective mass
419   //
420
421   Float_t mass = -0.1;
422   Bool_t sign = CheckSigns(v0);
423   if(0 == id){
424     mass = v0->GetEffMass(0, 0);
425   }
426   else if(1 == id){
427     mass = v0->GetEffMass(2, 2);
428   }
429   else if(2 == id){
430     mass = (sign) ? v0->GetEffMass(4, 2) : v0->GetEffMass(2, 4);
431   }
432   else if(3 == id){
433     mass = (sign) ? v0->GetEffMass(2, 4) : v0->GetEffMass(4, 2);
434   }
435   else{
436     AliWarning(Form("Unrecognized V0 id: %i", id));
437   }
438
439   return mass;
440
441 }
442 //__________________________________________________________
443 Bool_t AliAnalysisTaskCheckV0tender::CheckSigns(AliESDv0 * const v0){
444   //
445   // check wheter the sign was correctly applied to 
446   // V0 daughter tracks
447   // This function should become obsolete once the V0 finder will be updated (fixed)
448   //
449   
450   Bool_t correct = kFALSE;
451
452   Int_t pIndex = 0, nIndex = 0;
453   pIndex = v0->GetPindex();
454   nIndex = v0->GetNindex();
455   
456   AliESDtrack* d[2];
457   d[0] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(pIndex));
458   d[1] = dynamic_cast<AliESDtrack*>(fInputEvent->GetTrack(nIndex));
459
460   Int_t sign[2];
461   sign[0] = (int)d[0]->GetSign();
462   sign[1] = (int)d[1]->GetSign();
463   
464   if(-1 == sign[0] && 1 == sign[1]){
465     correct = kFALSE;
466     //v0->SetIndex(0, pIndex);  // set the index of the negative v0 track
467     //v0->SetIndex(1, nIndex);  // set the index of the positive v0 track
468   }
469   else{
470     correct = kTRUE;
471   }
472   
473   //pIndex = v0->GetPindex();
474   //nIndex = v0->GetNindex();
475   //printf("-D2: P: %i, N: %i\n", pIndex, nIndex);
476
477   return correct;
478 }