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