Adding task for like-sign bkg for D0->Kpi (Carmelo)
[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 fTotPosPairs(0),
58 fTotNegPairs(0),
59 fLsNormalization(1.)
60 {
61   //
62   // Default constructor
63   //
64 }
65
66 //________________________________________________________________________
67 AliAnalysisTaskSEBkgLikeSignD0::AliAnalysisTaskSEBkgLikeSignD0(const char *name):
68 AliAnalysisTaskSE(name),
69 fOutput(0),
70 fHistMassD0(0),
71 fHistMassLS(0),
72 fHistCtsD0(0),
73 fHistCtsLS(0),
74 fHistCtsLSpos(0),
75 fHistCtsLSneg(0),
76 fHistCPtaD0(0),
77 fHistCPtaLS(0),
78 fHistd0d0D0(0),
79 fHistd0d0LS(0),
80 fHistDCAD0(0),
81 fHistDCALS(0),
82 fVHF(0),
83 fTotPosPairs(0),
84 fTotNegPairs(0),
85 fLsNormalization(1.)
86 {
87   //
88   // Standard constructor
89   //
90   // Output slot #1 writes into a TList container
91   DefineOutput(1,TList::Class());  //My private output
92 }
93
94 //________________________________________________________________________
95 AliAnalysisTaskSEBkgLikeSignD0::~AliAnalysisTaskSEBkgLikeSignD0()
96 {
97   // Destructor
98   if (fOutput) {
99     delete fOutput;
100     fOutput = 0;
101   }
102
103   if (fVHF) {
104     delete fVHF;
105     fVHF = 0;
106   }
107
108 }
109 //________________________________________________________________________
110 void AliAnalysisTaskSEBkgLikeSignD0::Init()
111 {
112   // Initialization
113
114   if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignD0::Init() \n");
115
116   gROOT->LoadMacro("ConfigVertexingHF.C");
117
118   fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
119   fVHF->PrintStatus();
120
121   return;
122 }
123
124 //________________________________________________________________________
125 void AliAnalysisTaskSEBkgLikeSignD0::UserCreateOutputObjects()
126 {
127   // Create the output container
128   //
129   if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignD0::UserCreateOutputObjects() \n");
130
131   // Several histograms are more conveniently managed in a TList
132   fOutput = new TList();
133   fOutput->SetOwner();
134
135   fHistMassD0 = new TH1F("fHistMassD0", "D0 invariant mass; M [GeV]; Entries",200,1.56,2.16);
136   fHistMassD0->Sumw2();
137   fHistMassD0->SetMinimum(0);
138   fOutput->Add(fHistMassD0);
139
140   fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; M [GeV]; Entries",200,1.56,2.16);
141   fHistMassLS->Sumw2();
142   fHistMassLS->SetMinimum(0);
143   fOutput->Add(fHistMassLS);
144
145   fHistCtsD0 = new TH1F("fHistCtsD0", "D0 cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
146   fHistCtsD0->Sumw2();
147   fHistCtsD0->SetMinimum(0);
148   fOutput->Add(fHistCtsD0);
149
150   fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
151   fHistCtsLS->Sumw2();
152   fHistCtsLS->SetMinimum(0);
153   fOutput->Add(fHistCtsLS);
154
155   fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
156   fHistCtsLSpos->Sumw2();
157   fHistCtsLSpos->SetMinimum(0);
158   fOutput->Add(fHistCtsLSpos);
159
160   fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
161   fHistCtsLSneg->Sumw2();
162   fHistCtsLSneg->SetMinimum(0);
163   fOutput->Add(fHistCtsLSneg);
164
165   fHistCPtaD0 = new TH1F("fHistCPtaD0", "D0 cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.);
166   fHistCPtaD0->Sumw2();
167   fHistCPtaD0->SetMinimum(0);
168   fOutput->Add(fHistCPtaD0);
169
170   fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.);
171   fHistCPtaLS->Sumw2();
172   fHistCPtaLS->SetMinimum(0);
173   fOutput->Add(fHistCPtaLS);
174
175   fHistd0d0D0 = new TH1F("fHistd0d0D0", "D0 product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.);
176   fHistd0d0D0->Sumw2(); 
177   fHistd0d0D0->SetMinimum(0);
178   fOutput->Add(fHistd0d0D0);
179
180   fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.);
181   fHistd0d0LS->Sumw2();
182   fHistd0d0LS->SetMinimum(0);
183   fOutput->Add(fHistd0d0LS);
184
185   fHistDCAD0 = new TH1F("fHistDCAD0", "D0 distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
186   fHistDCAD0->Sumw2(); 
187   fHistDCAD0->SetMinimum(0);
188   fOutput->Add(fHistDCAD0);
189
190   fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
191   fHistDCALS->Sumw2(); 
192   fHistDCALS->SetMinimum(0);
193   fOutput->Add(fHistDCALS);
194
195   return;
196 }
197
198 //________________________________________________________________________
199 void AliAnalysisTaskSEBkgLikeSignD0::UserExec(Option_t */*option*/)
200 {
201   // Execute analysis for current event:
202   // heavy flavor candidates association to MC truth
203   
204   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
205
206   // load heavy flavour vertices
207   TClonesArray *arrayVerticesHF =
208     (TClonesArray*)aod->GetList()->FindObject("VerticesHF");
209   if(!arrayVerticesHF) {
210     printf("AliAnalysisTaskSEBkgLikeSignD0::UserExec: VerticesHF branch not found!\n");
211     return;
212   }
213
214   // load D0->ee candidates                                                   
215   TClonesArray *arrayD0toKpi =
216     (TClonesArray*)aod->GetList()->FindObject("D0toKpi");
217   if(!arrayD0toKpi) {
218     printf("AliAnalysisTaskSEBkgLikeSignD0::UserExec: D0toKpi branch not found!\n");
219     return;
220   }
221
222   // load like sign candidates
223   TClonesArray *arrayLikeSign =
224     (TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong");
225   if(!arrayLikeSign) {
226     printf("AliAnalysisTaskSEBkgLikeSignD0::UserExec: LikeSign2Prong branch not found!\n");
227     return;
228   }
229
230   // AOD primary vertex
231   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
232
233   // make trkIDtoEntry register (temporary)
234   Int_t trkIDtoEntry[100000];
235   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
236     AliAODTrack *track = aod->GetTrack(it);
237     trkIDtoEntry[track->GetID()]=it;
238   }
239
240   // loop over Like sign candidates
241   Int_t nPosPairs=0,nNegPairs=0;
242   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
243   printf("+++\n+++Number of like sign pairs ---> %d \n+++\n", nLikeSign);
244
245   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
246     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
247     Bool_t unsetvtx=kFALSE;
248     if(!d->GetOwnPrimaryVtx()) {
249         d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
250         unsetvtx=kTRUE;
251     }
252     Int_t okD0ls=0; Int_t okD0barls=0;
253     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0ls,okD0barls)) {
254        if(okD0ls)fHistMassLS->Fill(d->InvMassD0(),0.5);
255        if(okD0barls)fHistMassLS->Fill(d->InvMassD0bar(),0.5);
256        fHistCPtaLS->Fill(d->CosPointingAngle());
257        fHistd0d0LS->Fill(1e8*d->Prodd0d0());
258        if(okD0ls)fHistCtsLS->Fill(d->CosThetaStarD0(),0.5);
259        if(okD0barls)fHistCtsLS->Fill(d->CosThetaStarD0bar(),0.5);
260        fHistDCALS->Fill(100*d->GetDCA());
261        PostData(1,fOutput);
262        AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
263        if(!trk0) {
264           trk0=aod->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
265           printf("references to standard AOD not available \n");
266        }
267        if((trk0->Charge())==1) {
268           nPosPairs++;
269           fHistCtsLSpos->Fill(d->CosThetaStarD0());
270           PostData(1,fOutput);
271         } else {
272           nNegPairs++;
273           fHistCtsLSneg->Fill(d->CosThetaStarD0());
274           PostData(1,fOutput);
275         }
276       PostData(1,fOutput);
277       }
278     if(unsetvtx) d->UnsetOwnPrimaryVtx();
279   }
280
281   printf("------------ N. of positive pairs in Event ----- %d \n", nPosPairs);
282   printf("------------ N. of negative pairs in Event ----- %d \n", nNegPairs);
283
284   fTotPosPairs += nPosPairs;
285   fTotNegPairs += nNegPairs;
286
287   // loop over D0 candidates
288   Int_t nD0toKpi = arrayD0toKpi->GetEntriesFast();
289   printf("Number of like D0 -> Kpi candidates ---> %d \n", nD0toKpi);
290
291   for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
292     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->UncheckedAt(iD0toKpi);
293     Bool_t unsetvtx=kFALSE;
294     if(!d->GetOwnPrimaryVtx()) {
295       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
296       unsetvtx=kTRUE;
297     }
298     Int_t okD0=0; Int_t okD0bar=0;
299     if(d->SelectD0(fVHF->GetD0toKpiCuts(),okD0,okD0bar)) {
300       if(okD0)fHistMassD0->Fill(d->InvMassD0(),0.5);
301       if(okD0bar)fHistMassD0->Fill(d->InvMassD0bar(),0.5);
302       if(okD0)fHistCtsD0->Fill(d->CosThetaStarD0(),0.5);
303       if(okD0bar)fHistCtsD0->Fill(d->CosThetaStarD0bar(),0.5);
304       fHistd0d0D0->Fill(1e8*d->Prodd0d0());
305       fHistCPtaD0->Fill(d->CosPointingAngle());
306       fHistDCAD0->Fill(100*d->GetDCA());
307       PostData(1,fOutput);
308     }
309    if(unsetvtx) d->UnsetOwnPrimaryVtx();
310   }
311
312   return;
313 }
314
315 //________________________________________________________________________
316 void AliAnalysisTaskSEBkgLikeSignD0::Terminate(Option_t */*option*/)
317 {
318   // Terminate analysis
319   //
320   if(fDebug > 1) printf("AnalysisTaskSEBkgLikeSignD0: Terminate() \n");
321
322   fOutput = dynamic_cast<TList*> (GetOutputData(1));
323   if (!fOutput) {     
324     printf("ERROR: fOutput not available\n");
325     return;
326   }
327
328   fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs); 
329
330   fHistMassD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassD0"));
331   fHistMassLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
332   fHistCtsD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsD0"));
333   fHistCtsLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
334   fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
335   fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
336   fHistCPtaD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaD0"));
337   fHistCPtaLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
338   fHistd0d0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0D0"));
339   fHistd0d0LS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
340   fHistDCAD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCAD0"));
341   fHistDCALS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCALS"));
342
343   if(fLsNormalization>0) {
344     fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries());
345     fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries());
346     fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries());
347     fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries());
348     fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries());
349     fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries());
350     fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries());
351   }
352
353   return;
354 }