Use AliAnalysisTaskSE::AODEvent() to read the AOD that has just been produced by...
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEBkgLikeSignD0.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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 ///////////////////////////////////////////////////////////////////////////
17 //
18 // AliAnalysisTaskSE for reading both reconstructed D0->Kpi candidates
19 // and like sign pairs and for drawing corresponding distributions
20 //
21 // Author: C.Di Giglio, carmelo.digiglio@ba.infn.it
22 ///////////////////////////////////////////////////////////////////////////
23
24 #include <TSystem.h>
25 #include <TROOT.h>
26 #include <TClonesArray.h>
27 #include <TList.h>
28 #include <TH1F.h>
29
30 #include "AliAODEvent.h"
31 #include "AliAODVertex.h"
32 #include "AliAODTrack.h"
33 #include "AliAODRecoDecayHF2Prong.h"
34 #include "AliAnalysisVertexingHF.h"
35 #include "AliAnalysisTaskSE.h"
36 #include "AliAnalysisTaskSEBkgLikeSignD0.h"
37
38 ClassImp(AliAnalysisTaskSEBkgLikeSignD0)
39
40 //________________________________________________________________________
41 AliAnalysisTaskSEBkgLikeSignD0::AliAnalysisTaskSEBkgLikeSignD0():
42 AliAnalysisTaskSE(),
43 fOutput(0), 
44 fHistMassD0(0),
45 fHistMassLS(0),
46 fHistCtsD0(0),           
47 fHistCtsLS(0),
48 fHistCtsLSpos(0),
49 fHistCtsLSneg(0),
50 fHistCPtaD0(0),          
51 fHistCPtaLS(0),
52 fHistd0d0D0(0),          
53 fHistd0d0LS(0),
54 fHistDCAD0(0),           
55 fHistDCALS(0),
56 fVHF(0),
57 fNentries(0),
58 fTotPosPairs(0),
59 fTotNegPairs(0),
60 fLsNormalization(1.)
61 {
62   //
63   // Default constructor
64   //
65 }
66
67 //________________________________________________________________________
68 AliAnalysisTaskSEBkgLikeSignD0::AliAnalysisTaskSEBkgLikeSignD0(const char *name):
69 AliAnalysisTaskSE(name),
70 fOutput(0),
71 fHistMassD0(0),
72 fHistMassLS(0),
73 fHistCtsD0(0),
74 fHistCtsLS(0),
75 fHistCtsLSpos(0),
76 fHistCtsLSneg(0),
77 fHistCPtaD0(0),
78 fHistCPtaLS(0),
79 fHistd0d0D0(0),
80 fHistd0d0LS(0),
81 fHistDCAD0(0),
82 fHistDCALS(0),
83 fVHF(0),
84 fNentries(0),
85 fTotPosPairs(0),
86 fTotNegPairs(0),
87 fLsNormalization(1.)
88 {
89   //
90   // Standard constructor
91   //
92   // Output slot #1 writes into a TList container
93   DefineOutput(1,TList::Class());  //My private output
94   // Output slot #2 writes into a TH1F container
95   DefineOutput(2,TH1F::Class());  //My private output
96
97 }
98
99 //________________________________________________________________________
100 AliAnalysisTaskSEBkgLikeSignD0::~AliAnalysisTaskSEBkgLikeSignD0()
101 {
102   // Destructor
103   if (fOutput) {
104     delete fOutput;
105     fOutput = 0;
106   }
107
108   if (fVHF) {
109     delete fVHF;
110     fVHF = 0;
111   }
112
113   if (fNentries){
114     delete fNentries;
115     fNentries = 0;
116   }
117
118 }
119 //________________________________________________________________________
120 void AliAnalysisTaskSEBkgLikeSignD0::Init()
121 {
122   // Initialization
123
124   if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignD0::Init() \n");
125
126   gROOT->LoadMacro("ConfigVertexingHF.C");
127
128   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
129   fVHF->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,1,-0.00025,0.8);
130   fVHF->PrintStatus();
131
132   return;
133 }
134
135 //________________________________________________________________________
136 void AliAnalysisTaskSEBkgLikeSignD0::UserCreateOutputObjects()
137 {
138   // Create the output container
139   //
140   if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignD0::UserCreateOutputObjects() \n");
141
142   // Several histograms are more conveniently managed in a TList
143   fOutput = new TList();
144   fOutput->SetOwner();
145
146   fHistMassD0 = new TH1F("fHistMassD0", "D0 invariant mass; M [GeV]; Entries",200,1.765,1.965);
147   fHistMassD0->Sumw2();
148   fHistMassD0->SetMinimum(0);
149   fOutput->Add(fHistMassD0);
150
151   fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; M [GeV]; Entries",200,1.765,1.965);
152   fHistMassLS->Sumw2();
153   fHistMassLS->SetMinimum(0);
154   fOutput->Add(fHistMassLS);
155
156   fHistCtsD0 = new TH1F("fHistCtsD0", "D0 cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
157   fHistCtsD0->Sumw2();
158   fHistCtsD0->SetMinimum(0);
159   fOutput->Add(fHistCtsD0);
160
161   fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
162   fHistCtsLS->Sumw2();
163   fHistCtsLS->SetMinimum(0);
164   fOutput->Add(fHistCtsLS);
165
166   fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
167   fHistCtsLSpos->Sumw2();
168   fHistCtsLSpos->SetMinimum(0);
169   fOutput->Add(fHistCtsLSpos);
170
171   fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
172   fHistCtsLSneg->Sumw2();
173   fHistCtsLSneg->SetMinimum(0);
174   fOutput->Add(fHistCtsLSneg);
175
176   fHistCPtaD0 = new TH1F("fHistCPtaD0", "D0 cosine of pointing angle; Cos#Theta_{point}; Entries",200,0,1.);
177   fHistCPtaD0->Sumw2();
178   fHistCPtaD0->SetMinimum(0);
179   fOutput->Add(fHistCPtaD0);
180
181   fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,0,1.);
182   fHistCPtaLS->Sumw2();
183   fHistCPtaLS->SetMinimum(0);
184   fOutput->Add(fHistCPtaLS);
185
186   fHistd0d0D0 = new TH1F("fHistd0d0D0", "D0 product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.);
187   fHistd0d0D0->Sumw2(); 
188   fHistd0d0D0->SetMinimum(0);
189   fOutput->Add(fHistd0d0D0);
190
191   fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.);
192   fHistd0d0LS->Sumw2();
193   fHistd0d0LS->SetMinimum(0);
194   fOutput->Add(fHistd0d0LS);
195
196   fHistDCAD0 = new TH1F("fHistDCAD0", "D0 distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
197   fHistDCAD0->Sumw2(); 
198   fHistDCAD0->SetMinimum(0);
199   fOutput->Add(fHistDCAD0);
200
201   fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
202   fHistDCALS->Sumw2(); 
203   fHistDCALS->SetMinimum(0);
204   fOutput->Add(fHistDCALS);
205
206   fNentries=new TH1F("nentriesLS", "Look at the number of entries! it is = to the  number of AODs", 2,1.,2.);
207
208   return;
209 }
210
211 //________________________________________________________________________
212 void AliAnalysisTaskSEBkgLikeSignD0::UserExec(Option_t */*option*/)
213 {
214   // Execute analysis for current event:
215   // heavy flavor candidates association to MC truth
216   
217   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
218
219   // In case there is an AOD handler writing a standard AOD, use the AOD 
220   // event in memory rather than the input (ESD) event.
221   if (!aod && AODEvent() && IsStandardAOD()) aod = dynamic_cast<AliAODEvent*> (AODEvent());
222
223   // load heavy flavour vertices
224   TClonesArray *arrayVerticesHF =
225     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
226   if(!arrayVerticesHF) {
227     printf("AliAnalysisTaskSEBkgLikeSignD0::UserExec: VerticesHF branch not found!\n");
228     return;
229   }
230
231   // load D0->ee candidates                                                   
232   TClonesArray *arrayD0toKpi =
233     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
234   if(!arrayD0toKpi) {
235     printf("AliAnalysisTaskSEBkgLikeSignD0::UserExec: D0toKpi branch not found!\n");
236     return;
237   }
238
239   // load like sign candidates
240   TClonesArray *arrayLikeSign =
241     (TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong");
242   if(!arrayLikeSign) {
243     printf("AliAnalysisTaskSEBkgLikeSignD0::UserExec: LikeSign2Prong branch not found!\n");
244     return;
245   }
246
247   // AOD primary vertex
248   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
249
250   // make trkIDtoEntry register (temporary)
251   Int_t trkIDtoEntry[100000];
252   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
253     AliAODTrack *track = aod->GetTrack(it);
254     trkIDtoEntry[track->GetID()]=it;
255   }
256
257   //histogram filled with 1 for every AOD
258   fNentries->Fill(1);
259   PostData(2,fNentries);
260
261   // loop over Like sign candidates
262   Int_t nPosPairs=0,nNegPairs=0;
263   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
264   if(fDebug>1) printf("+++\n+++Number of like sign pairs ---> %d \n+++\n", nLikeSign);
265
266   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
267     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
268     Bool_t unsetvtx=kFALSE;
269     if(!d->GetOwnPrimaryVtx()) {
270         d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
271         unsetvtx=kTRUE;
272     }
273     Int_t okD0ls=0; Int_t okD0barls=0;
274     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0ls,okD0barls)) {
275        AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
276        AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
277        if(!trk0 || !trk1) {
278           trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
279           trk1=aod->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
280           printf("references to standard AOD not available \n");
281        }
282        if(okD0ls) fHistMassLS->Fill(d->InvMassD0());
283        if(okD0barls) fHistMassLS->Fill(d->InvMassD0bar());
284        fHistCPtaLS->Fill(d->CosPointingAngle());
285        fHistd0d0LS->Fill(1e8*d->Prodd0d0());
286        if(okD0ls) fHistCtsLS->Fill(d->CosThetaStarD0());
287        if(okD0barls) fHistCtsLS->Fill(d->CosThetaStarD0bar());
288        fHistDCALS->Fill(100*d->GetDCA());
289        //PostData(1,fOutput);
290        if((trk0->Charge())==1) {
291           nPosPairs++;
292           fHistCtsLSpos->Fill(d->CosThetaStarD0());
293           //PostData(1,fOutput);
294         } else {
295           nNegPairs++;
296           fHistCtsLSneg->Fill(d->CosThetaStarD0());
297           //PostData(1,fOutput);
298         }
299       PostData(1,fOutput);
300       }
301     if(unsetvtx) d->UnsetOwnPrimaryVtx();
302   }
303
304   if(fDebug>1) printf("------------ N. of positive pairs in Event ----- %d \n", nPosPairs);
305   if(fDebug>1) printf("------------ N. of negative pairs in Event ----- %d \n", nNegPairs);
306
307   fTotPosPairs += nPosPairs;
308   fTotNegPairs += nNegPairs;
309
310   // loop over D0 candidates
311   Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
312   if(fDebug>1) printf("Number of like D0 -> Kpi candidates ---> %d \n", nD0toKpi);
313
314   for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
315     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
316     Bool_t unsetvtx=kFALSE;
317     if(!d->GetOwnPrimaryVtx()) {
318       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
319       unsetvtx=kTRUE;
320     }
321     Int_t okD0=0; Int_t okD0bar=0;
322     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
323       fHistMassD0->Fill(d->InvMassD0());
324       fHistMassD0->Fill(d->InvMassD0bar());
325       fHistCtsD0->Fill(d->CosThetaStarD0());
326       fHistCtsD0->Fill(d->CosThetaStarD0bar());
327       fHistd0d0D0->Fill(1e8*d->Prodd0d0());
328       fHistCPtaD0->Fill(d->CosPointingAngle());
329       fHistDCAD0->Fill(100*d->GetDCA());
330       PostData(1,fOutput);
331     }
332    if(unsetvtx) d->UnsetOwnPrimaryVtx();
333   }
334
335   return;
336 }
337
338 //________________________________________________________________________
339 void AliAnalysisTaskSEBkgLikeSignD0::Terminate(Option_t */*option*/)
340 {
341   // Terminate analysis
342   //
343   if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignD0: Terminate() \n");
344
345   fOutput = dynamic_cast<TList*> (GetOutputData(1));
346   if (!fOutput) {     
347     printf("ERROR: fOutput not available\n");
348     return;
349   }
350
351   fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs); 
352
353   fHistMassD0   = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassD0"));
354   fHistMassLS   = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
355   fHistCtsD0    = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsD0"));
356   fHistCtsLS    = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
357   fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
358   fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
359   fHistCPtaD0   = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaD0"));
360   fHistCPtaLS   = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
361   fHistd0d0D0   = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0D0"));
362   fHistd0d0LS   = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
363   fHistDCAD0    = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCAD0"));
364   fHistDCALS    = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCALS"));
365
366   if(fLsNormalization>0) {
367     fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries());
368     fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries());
369     fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries());
370     fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries());
371     fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries());
372     fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries());
373     fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries());
374   }
375
376   return;
377 }