Updates in PbPb cuts (Andrea Rossi)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSECompareHF.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, 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 //
20 // AliAnalysisTaskSE for the comparison of heavy flavor
21 // decay candidates with the MC truth.
22 //
23 // Author: A.Dainese, andrea.dainese@lnl.infn.it
24 /////////////////////////////////////////////////////////////
25
26 #include <TClonesArray.h>
27 #include <TNtuple.h>
28 #include <TList.h>
29 #include <TH1F.h>
30
31 #include "AliAnalysisManager.h"
32 #include "AliAODHandler.h"
33 #include "AliAODEvent.h"
34 #include "AliAODVertex.h"
35 #include "AliESDtrack.h"
36 #include "AliAODTrack.h"
37 #include "AliAODMCHeader.h"
38 #include "AliAODMCParticle.h"
39 #include "AliAODRecoDecayHF2Prong.h"
40 #include "AliAODRecoDecayHF3Prong.h"
41 #include "AliAODRecoDecayHF4Prong.h"
42 #include "AliAODRecoCascadeHF.h"
43 #include "AliAnalysisVertexingHF.h"
44 #include "AliAnalysisTaskSE.h"
45 #include "AliAnalysisTaskSECompareHF.h"
46
47 ClassImp(AliAnalysisTaskSECompareHF)
48
49
50 //________________________________________________________________________
51 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF():
52 AliAnalysisTaskSE(),
53 fOutput(0), 
54 fNtupleCmp(0),
55 fHistMass(0),
56 fHistNEvents(0)
57 {
58   // Default constructor
59
60   // NO DefineOutput() HERE (ONLY IN STANDARD CONSTRUCTOR)
61 }
62
63 //________________________________________________________________________
64 AliAnalysisTaskSECompareHF::AliAnalysisTaskSECompareHF(const char *name):
65 AliAnalysisTaskSE(name),
66 fOutput(0), 
67 fNtupleCmp(0),
68 fHistMass(0),
69 fHistNEvents(0)
70 {
71   // Standard constructor
72
73   // Output slot #1 writes into a TList container
74   DefineOutput(1,TList::Class());  //My private output
75   // Output slot #2 writes into a TNtuple container
76   DefineOutput(2,TNtuple::Class());  //My private output
77 }
78
79 //________________________________________________________________________
80 AliAnalysisTaskSECompareHF::~AliAnalysisTaskSECompareHF()
81 {
82   // Destructor
83   if (fOutput) {
84     delete fOutput;
85     fOutput = 0;
86   }
87 }  
88
89 //________________________________________________________________________
90 void AliAnalysisTaskSECompareHF::Init()
91 {
92   // Initialization
93
94   if(fDebug > 1) printf("AnalysisTaskSECompareHF::Init() \n");
95   
96   return;
97 }
98
99 //________________________________________________________________________
100 void AliAnalysisTaskSECompareHF::UserCreateOutputObjects()
101 {
102   // Create the output container
103   //
104   if(fDebug > 1) printf("AnalysisTaskSECompareHF::UserCreateOutputObjects() \n");
105
106   // Several histograms are more conveniently managed in a TList
107   fOutput = new TList();
108   fOutput->SetOwner();
109
110   fHistMass = new TH1F("fHistMass", "D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
111   fHistMass->Sumw2();
112   fHistMass->SetMinimum(0);
113   fOutput->Add(fHistMass);
114
115   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events; ; Events",3,-1.5,1.5);
116   fHistNEvents->Sumw2();
117   fHistNEvents->SetMinimum(0);
118   fOutput->Add(fHistNEvents);
119
120   OpenFile(2); // 2 is the slot number of the ntuple
121   fNtupleCmp = new TNtuple("fNtupleCmp","Charm comparison","pdg:nprongs:VxRec:VxTrue:ErrVx:VyRec:VyTrue:ErrVy:VzRec:VzTrue:ErrVz:Chi2toNDF:PtRec:Mrec:CPta:Prodd0");
122
123   return;
124 }
125
126 //________________________________________________________________________
127 void AliAnalysisTaskSECompareHF::UserExec(Option_t */*option*/)
128 {
129   // Execute analysis for current event:
130   // heavy flavor candidates association to MC truth
131
132   
133   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
134
135   TClonesArray *inputArrayVertices = 0;
136   TClonesArray *inputArrayD0toKpi = 0;
137   TClonesArray *inputArrayDstar = 0;
138
139   if(!aod && AODEvent() && IsStandardAOD()) {
140     // In case there is an AOD handler writing a standard AOD, use the AOD 
141     // event in memory rather than the input (ESD) event.    
142     aod = dynamic_cast<AliAODEvent*> (AODEvent());
143     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
144     // have to taken from the AOD event hold by the AliAODExtension
145     AliAODHandler* aodHandler = (AliAODHandler*) 
146       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
147     if(aodHandler->GetExtensions()) {
148       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
149       AliAODEvent *aodFromExt = ext->GetAOD();
150       // load HF vertices                
151       inputArrayVertices = (TClonesArray*)aodFromExt->GetList()->FindObject("VerticesHF");
152       // load D0->Kpi candidates
153       inputArrayD0toKpi = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
154       // load D*+ candidates                                                   
155       inputArrayDstar = (TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
156     }
157   } else if(aod) {
158     // load HF vertices                
159     inputArrayVertices = (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
160     // load D0->Kpi candidates                                                 
161     inputArrayD0toKpi = (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
162     // load D*+ candidates                                                   
163     inputArrayDstar = (TClonesArray*)aod->GetList()->FindObject("Dstar");
164   }
165
166
167   if(!inputArrayVertices || !aod) {
168     printf("AliAnalysisTaskSECompareHF::UserExec: Vertices branch not found!\n");
169     return;
170   }
171   if(!inputArrayD0toKpi) {
172     printf("AliAnalysisTaskSECompareHF::UserExec: D0toKpi branch not found!\n");
173     return;
174   }
175   if(!inputArrayDstar) {
176     printf("AliAnalysisTaskSECompareHF::UserExec: Dstar branch not found!\n");
177     return;
178   }
179   
180
181   fHistNEvents->Fill(0); // count event
182   // Post the data already here
183   PostData(1,fOutput);
184
185   // AOD primary vertex
186   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
187   //vtx1->Print();
188
189   // load MC particles
190   TClonesArray *mcArray = 
191     (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
192   if(!mcArray) {
193     printf("AliAnalysisTaskSECompareHF::UserExec: MC particles branch not found!\n");
194     return;
195   }
196
197   // load MC header
198   AliAODMCHeader *mcHeader = 
199     (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
200   if(!mcHeader) {
201     printf("AliAnalysisTaskSECompareHF::UserExec: MC header branch not found!\n");
202     return;
203   }
204   
205   Int_t nprongs,lab,okD0,okD0bar,pdg;
206   Bool_t unsetvtx;    
207   Double_t invmass,posRec[3],posTrue[3],covRec[6],errx,erry,errz;
208   AliAODRecoDecayHF2Prong *d2=0;
209   AliAODRecoDecayHF3Prong *d3=0;
210   AliAODRecoDecayHF4Prong *d4=0;
211
212   Int_t pdgDgD0toKpi[2]={321,211};
213   Int_t pdgDgDplustoKpipi[3]={321,211,211};
214   Int_t pdgDgD0toKpipipi[4]={321,211,211,211};
215
216   // loop over vertices
217   Int_t nVertices = inputArrayVertices->GetEntriesFast();
218   if(fDebug>1) printf("Number of vertices: %d\n",nVertices);
219
220   for (Int_t iVtx = 0; iVtx < nVertices; iVtx++) {
221     AliAODVertex *vtx = (AliAODVertex*)inputArrayVertices->UncheckedAt(iVtx);
222
223     vtx->GetXYZ(posRec);  
224     vtx->GetCovarianceMatrix(covRec);
225     errx=1.; erry=1.; errz=1.;
226     if(covRec[0]>0) errx = TMath::Sqrt(covRec[0]);
227     if(covRec[2]>0) erry = TMath::Sqrt(covRec[2]);
228     if(covRec[5]>0) errz = TMath::Sqrt(covRec[5]);
229           
230
231     // get parent AliAODRecoDecayHF
232     nprongs= vtx->GetNDaughters();
233     // check that the daughters have kITSrefit and kTPCrefit
234     Bool_t allDgOK=kTRUE;
235     for(Int_t idg=0; idg<nprongs; idg++) {
236       AliAODTrack *track = (AliAODTrack*)vtx->GetDaughter(idg);
237       if(!(track->GetStatus()&AliESDtrack::kITSrefit)) allDgOK=kFALSE;
238       if(!(track->GetStatus()&AliESDtrack::kTPCrefit)) allDgOK=kFALSE;
239     }
240     if(!allDgOK) continue;
241
242
243     switch(nprongs) {
244     case 2: // look for D0->Kpi
245       d2 = (AliAODRecoDecayHF2Prong*)vtx->GetParent();
246       if(d2->IsLikeSign()) continue;
247       if(d2->Charge() != 0) continue; // these are D* 
248       lab = d2->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
249       if(lab>=0) {
250         unsetvtx=kFALSE;
251         if(!d2->GetOwnPrimaryVtx()) {
252           d2->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
253           unsetvtx=kTRUE;
254         }
255         okD0=1; okD0bar=1; 
256         AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
257         pdg = dMC->GetPdgCode();
258         invmass = (pdg==421 ? d2->InvMassD0() : d2->InvMassD0bar());
259         // get a daughter for true pos of decay vertex
260         AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
261         dg0MC->XvYvZv(posTrue);
262         fHistMass->Fill(invmass);
263         // Post the data already here
264         PostData(1,fOutput);
265         Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
266                          (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
267                          (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
268                          (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
269                          (Float_t)d2->GetReducedChi2(),(Float_t)d2->Pt(),(Float_t)invmass,
270                          (Float_t)d2->CosPointingAngle(),(Float_t)d2->Prodd0d0()};
271         fNtupleCmp->Fill(tmp);
272         PostData(2,fNtupleCmp);
273       
274         if(unsetvtx) d2->UnsetOwnPrimaryVtx();
275       }
276       break;
277     case 3: // look for D+
278       d3 = (AliAODRecoDecayHF3Prong*)vtx->GetParent();
279       if(d3->IsLikeSign()) continue;
280       lab = d3->MatchToMC(411,mcArray,3,pdgDgDplustoKpipi);
281       if(lab>=0) {
282         unsetvtx=kFALSE;
283         if(!d3->GetOwnPrimaryVtx()) {
284           d3->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
285           unsetvtx=kTRUE;
286         }
287         AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
288         pdg = dMC->GetPdgCode();
289         invmass = d3->InvMassDplus();
290         // get a daughter for true pos of decay vertex
291         AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
292         dg0MC->XvYvZv(posTrue);
293         Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
294                          (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
295                          (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
296                          (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
297                          (Float_t)d3->GetReducedChi2(),(Float_t)d3->Pt(),(Float_t)invmass,
298                          (Float_t)d3->CosPointingAngle(),(Float_t)(d3->Getd0Prong(0)*d3->Getd0Prong(1)*d3->Getd0Prong(2))};
299         fNtupleCmp->Fill(tmp);
300         PostData(2,fNtupleCmp);
301         if(unsetvtx) d3->UnsetOwnPrimaryVtx();
302       }
303       break;
304     case 4: // look for D0->Kpipipi
305       d4 = (AliAODRecoDecayHF4Prong*)vtx->GetParent();
306       if(d4->IsLikeSign()) continue;
307       lab = d4->MatchToMC(421,mcArray,4,pdgDgD0toKpipipi);
308       if(lab>=0) {
309         unsetvtx=kFALSE;
310         if(!d4->GetOwnPrimaryVtx()) {
311           d4->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
312           unsetvtx=kTRUE;
313         }
314         okD0=0; okD0bar=0; 
315         AliAODMCParticle *dMC = (AliAODMCParticle*)mcArray->At(lab);
316         pdg = dMC->GetPdgCode();
317         //invmass = (pdg==421 ? d->InvMassD0() : d->InvMassD0bar());
318         invmass = 10.;    // dummy
319         // get a daughter for true pos of decay vertex
320         AliAODMCParticle *dg0MC = (AliAODMCParticle*)mcArray->At(dMC->GetDaughter(0));
321         dg0MC->XvYvZv(posTrue);
322         Float_t tmp[16]={(Float_t)pdg,(Float_t)nprongs,
323                          (Float_t)posRec[0],(Float_t)posTrue[0],(Float_t)errx,
324                          (Float_t)posRec[1],(Float_t)posTrue[1],(Float_t)erry,
325                          (Float_t)posRec[2],(Float_t)posTrue[2],(Float_t)errz,
326                          (Float_t)d4->GetReducedChi2(),(Float_t)d4->Pt(),(Float_t)invmass,
327                          (Float_t)d4->CosPointingAngle(),(Float_t)(d4->Getd0Prong(0)*d4->Getd0Prong(1)*d4->Getd0Prong(2)*d4->Getd0Prong(3))};
328         fNtupleCmp->Fill(tmp);
329         PostData(2,fNtupleCmp);
330         if(unsetvtx) d4->UnsetOwnPrimaryVtx();
331       }
332       break;
333     }
334
335   } // end loop on vertices
336
337
338   // loop over D*+ candidates
339   /*
340   for (Int_t iDstar = 0; iDstar < inputArrayDstar->GetEntries(); iDstar++) {
341     AliAODRecoCascadeHF *c = (AliAODRecoCascadeHF*)inputArrayDstar->UncheckedAt(iDstar);
342     Int_t labDstar = c->MatchToMC(413,421,mcArray);
343     if(labDstar>=0) printf("GOOD MATCH FOR D*+\n"); 
344   }
345   */
346
347   // Post the data already here
348   PostData(1,fOutput);
349   PostData(2,fNtupleCmp);
350
351   return;
352 }
353 //________________________________________________________________________
354 void AliAnalysisTaskSECompareHF::Terminate(Option_t */*option*/)
355 {
356   // Terminate analysis
357   //
358   if(fDebug > 1) printf("AnalysisTaskSECompareHF: Terminate() \n");
359
360   fOutput = dynamic_cast<TList*> (GetOutputData(1));
361   if (!fOutput) {     
362     printf("ERROR: fOutput not available\n");
363     return;
364   }
365
366   fHistMass = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMass"));
367   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
368
369   //fNtupleCmp = dynamic_cast<TNtuple*> (GetOutputData(2));
370
371   return;
372 }
373