New class for normalization studies (Giacomo)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSED0Mass.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 D0 candidates invariant mass histogram
19 // and comparison with the MC truth and cut variables distributions.
20 //
21 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
22 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
23 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
24 /////////////////////////////////////////////////////////////
25
26 #include <Riostream.h>
27 #include <TClonesArray.h>
28 #include <TCanvas.h>
29 #include <TNtuple.h>
30 #include <TList.h>
31 #include <TH1F.h>
32 #include <TH2F.h>
33 #include <TDatabasePDG.h>
34
35 #include <AliAnalysisDataSlot.h>
36 #include <AliAnalysisDataContainer.h>
37 #include "AliAnalysisManager.h"
38 #include "AliESDtrack.h"
39 #include "AliVertexerTracks.h"
40 #include "AliAODHandler.h"
41 #include "AliAODEvent.h"
42 #include "AliAODVertex.h"
43 #include "AliAODTrack.h"
44 #include "AliAODMCHeader.h"
45 #include "AliAODMCParticle.h"
46 #include "AliAODRecoDecayHF2Prong.h"
47 #include "AliAODRecoCascadeHF.h"
48 #include "AliAnalysisVertexingHF.h"
49 #include "AliAnalysisTaskSE.h"
50 #include "AliAnalysisTaskSED0Mass.h"
51 #include "AliNormalizationCounter.h"
52
53 ClassImp(AliAnalysisTaskSED0Mass)
54
55
56 //________________________________________________________________________
57 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
58 AliAnalysisTaskSE(),
59 fOutputMass(0),
60 fDistr(0),
61 fNentries(0), 
62 fChecks(0),
63 fCuts(0),
64 fArray(0),
65 fReadMC(0),
66 fCutOnDistr(0),
67 fUsePid4Distr(0),
68 fCounter(0),
69 fNPtBins(1),
70 fTotPosPairs(0),
71 fTotNegPairs(0),
72 fLsNormalization(1.),
73 fFillOnlyD0D0bar(0)
74
75
76 {
77   // Default constructor
78 }
79
80 //________________________________________________________________________
81 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
82 AliAnalysisTaskSE(name),
83 fOutputMass(0), 
84 fDistr(0),
85 fNentries(0),
86 fChecks(0),
87 fCuts(0),
88 fArray(0),
89 fReadMC(0),
90 fCutOnDistr(0),
91 fUsePid4Distr(0),
92 fCounter(0),
93 fNPtBins(1),
94 fTotPosPairs(0),
95 fTotNegPairs(0),
96 fLsNormalization(1.),
97 fFillOnlyD0D0bar(0)
98
99 {
100   // Default constructor
101
102   fNPtBins=cuts->GetNPtBins();
103   fTotPosPairs=new Int_t[fNPtBins];
104   fTotNegPairs=new Int_t[fNPtBins];
105   for(Int_t i=0;i<fNPtBins;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
106     
107   fCuts=cuts;
108
109   // Output slot #1 writes into a TList container (mass with cuts)
110   DefineOutput(1,TList::Class());  //My private output
111   // Output slot #2 writes into a TList container (distributions)
112   DefineOutput(2,TList::Class());  //My private output
113   // Output slot #3 writes into a TH1F container (number of events)
114   DefineOutput(3,TH1F::Class());  //My private output
115   // Output slot #4 writes into a TList container (quality check)
116   DefineOutput(4,TList::Class());  //My private output
117   // Output slot #5 writes into a TList container (cuts)
118   DefineOutput(5,AliRDHFCutsD0toKpi::Class());  //My private output
119   // Output slot #6 writes Normalization Counter 
120   DefineOutput(6,AliNormalizationCounter::Class());
121 }
122
123 //________________________________________________________________________
124 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
125 {
126    if (fOutputMass) {
127     delete fOutputMass;
128     fOutputMass = 0;
129   }
130   if (fDistr) {
131     delete fDistr;
132     fDistr = 0;
133   }
134   if (fChecks) {
135     delete fChecks;
136     fChecks = 0;
137   }
138   if (fCuts) {
139     delete fCuts;
140     fCuts = 0;
141   }
142   if (fNentries){
143     delete fNentries;
144     fNentries = 0;
145   }
146   if(fCounter){
147     delete fCounter;
148     fCounter=0;
149   }
150  
151 }  
152
153 //________________________________________________________________________
154 void AliAnalysisTaskSED0Mass::Init()
155 {
156   // Initialization
157
158   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
159
160   
161   AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
162   const char* nameoutput=GetOutputSlot(5)->GetContainer()->GetName();
163   copyfCuts->SetName(nameoutput);
164   // Post the data
165   PostData(5,copyfCuts);
166
167
168   return;
169 }
170
171 //________________________________________________________________________
172 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
173 {
174
175   // Create the output container
176   //
177   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
178
179   // Several histograms are more conveniently managed in a TList
180   fOutputMass = new TList();
181   fOutputMass->SetOwner();
182   fOutputMass->SetName("listMass");
183
184   fDistr = new TList();
185   fDistr->SetOwner();
186   fDistr->SetName("distributionslist");
187
188   fChecks = new TList();
189   fChecks->SetOwner();
190   fChecks->SetName("checkHistograms");
191
192   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
193
194   for(Int_t i=0;i<fCuts->GetNPtBins();i++){
195
196     nameMass="histMass_";
197     nameMass+=i;
198     nameSgn27="histSgn27_";
199     nameSgn27+=i;
200     nameSgn="histSgn_";
201     nameSgn+=i;
202     nameBkg="histBkg_";
203     nameBkg+=i;
204     nameRfl="histRfl_";
205     nameRfl+=i;
206     nameMassNocutsS="hMassS_";
207     nameMassNocutsS+=i;
208     nameMassNocutsB="hMassB_";
209     nameMassNocutsB+=i;
210
211     //histograms of cut variable distributions
212
213     //  pT
214     namedistr="hptpiS_";
215     namedistr+=i;
216     TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
217
218     namedistr="hptKS_";
219     namedistr+=i;
220     TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
221
222     namedistr="hptB_";
223     namedistr+=i;
224     TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
225
226     //  pT no mass cut
227     namedistr="hptpiSnoMcut_";
228     namedistr+=i;
229     TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
230
231     namedistr="hptKSnoMcut_";
232     namedistr+=i;
233     TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
234
235     namedistr="hptB1prongnoMcut_";
236     namedistr+=i;
237     TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
238
239     namedistr="hptB2prongsnoMcut_";
240     namedistr+=i;
241     TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
242
243     
244     //  dca
245     namedistr="hdcaS_";
246     namedistr+=i;
247     TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
248     namedistr="hdcaB_";
249     namedistr+=i;
250     TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
251
252     //  costhetastar
253     namedistr="hcosthetastarS_";
254     namedistr+=i;
255     TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
256     namedistr="hcosthetastarB_";
257     namedistr+=i;
258     TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
259
260     // impact parameter
261     namedistr="hd0piS_";
262     namedistr+=i;
263     TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
264
265     namedistr="hd0KS_";
266     namedistr+=i;
267     TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
268
269     namedistr="hd0B_";
270     namedistr+=i;
271     TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
272
273     namedistr="hd0p0B_";
274     namedistr+=i;
275     TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
276
277     namedistr="hd0moresB_";
278     namedistr+=i;
279     TH1F *hd0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
280
281     namedistr="hd0d0moresB_";
282     namedistr+=i;
283     TH1F *hd0d0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
284
285     namedistr="hd0vmoresB_";
286     namedistr+=i;
287     TH1F *hd0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
288
289     namedistr="hd0d0vmoresB_";
290     namedistr+=i;
291     TH1F *hd0d0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
292
293     namedistr="hd0p1B_";
294     namedistr+=i;
295     TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
296
297     namedistr="hd0vpiS_";
298     namedistr+=i;
299     TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
300
301     namedistr="hd0vKS_";
302     namedistr+=i;
303     TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
304     namedistr="hd0vB_";
305     namedistr+=i;
306     TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
307
308     namedistr="hd0vp0B_";
309     namedistr+=i;
310     TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
311
312     namedistr="hd0vp1B_";
313     namedistr+=i;
314     TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
315
316     namedistr="hd0d0S_";
317     namedistr+=i;
318     TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
319     namedistr="hd0d0B_";
320     namedistr+=i;
321     TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
322
323      namedistr="hd0d0vS_";
324     namedistr+=i;
325     TH1F *hd0d0vS = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
326     namedistr="hd0d0vB_";
327     namedistr+=i;
328     TH1F *hd0d0vB = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
329
330     //decay lenght
331     namedistr="hdeclS_";
332     namedistr+=i;
333     TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length distribution;Decay Length [cm]",200,0,0.15);
334
335     namedistr="hdeclB_";
336     namedistr+=i;
337     TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length distribution;Decay Length [cm]",200,0,0.15);
338
339     namedistr="hnormdeclS_";
340     namedistr+=i;
341     TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;Decay Length/Err ",200,0,6.);
342
343     namedistr="hnormdeclB_";
344     namedistr+=i;
345     TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;Decay Length/Err ",200,0,6.);
346
347     namedistr="hdeclvS_";
348     namedistr+=i;
349     TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
350
351     namedistr="hdeclvB_";
352     namedistr+=i;
353     TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
354
355     namedistr="hnormdeclvS_";
356     namedistr+=i;
357     TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
358
359     namedistr="hnormdeclvB_";
360     namedistr+=i;
361     TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
362
363    //  costhetapoint
364     namedistr="hcosthetapointS_";
365     namedistr+=i;
366     TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
367     namedistr="hcosthetapointB_";
368     namedistr+=i;
369     TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
370
371     namedistr="hcosthetapointmoresB_";
372     namedistr+=i;
373     TH1F *hcosthetapointmoresB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
374
375    //  costhetapoint vs d0 or d0d0
376     namedistr="hcosthpointd0S_";
377     namedistr+=i;
378     TH2F *hcosthpointd0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
379
380     namedistr="hcosthpointd0B_";
381     namedistr+=i;
382     TH2F *hcosthpointd0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
383
384     namedistr="hcosthpointd0d0S_";
385     namedistr+=i;
386     TH2F *hcosthpointd0d0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
387     namedistr="hcosthpointd0d0B_";
388     namedistr+=i;
389     TH2F *hcosthpointd0d0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
390
391     fDistr->Add(hptpiS);
392     fDistr->Add(hptKS);
393     fDistr->Add(hptB);
394
395     fDistr->Add(hptpiSnoMcut);
396     fDistr->Add(hptKSnoMcut);
397     fDistr->Add(hptB1pnoMcut);
398     fDistr->Add(hptB2pnoMcut);
399
400     fDistr->Add(hdcaS);
401     fDistr->Add(hdcaB);
402
403     fDistr->Add(hd0piS);
404     fDistr->Add(hd0KS);
405     fDistr->Add(hd0B);
406     fDistr->Add(hd0p0B);
407     fDistr->Add(hd0moresB);
408     fDistr->Add(hd0p1B);
409
410     fDistr->Add(hd0vpiS);
411     fDistr->Add(hd0vKS);
412     fDistr->Add(hd0vB);
413     fDistr->Add(hd0vp0B);
414     fDistr->Add(hd0vp1B);
415     fDistr->Add(hd0vmoresB);
416
417     fDistr->Add(hd0d0S);
418     fDistr->Add(hd0d0B);
419     fDistr->Add(hd0d0moresB);
420
421     fDistr->Add(hd0d0vS);
422     fDistr->Add(hd0d0vB);
423     fDistr->Add(hd0d0vmoresB);
424
425     fDistr->Add(hcosthetastarS);
426     fDistr->Add(hcosthetastarB);
427
428     fDistr->Add(hcosthetapointS);
429     fDistr->Add(hcosthetapointB);
430     fDistr->Add(hcosthetapointmoresB);
431
432     fDistr->Add(hdeclengthS);
433     fDistr->Add(hdeclengthB);
434
435     fDistr->Add(hnormdeclengthS);
436     fDistr->Add(hnormdeclengthB);
437
438     fDistr->Add(hdeclengthvS);
439     fDistr->Add(hdeclengthvB);
440
441     fDistr->Add(hnormdeclengthvS);
442     fDistr->Add(hnormdeclengthvB);
443
444     fDistr->Add(hcosthpointd0S);
445     fDistr->Add(hcosthpointd0B);
446
447     fDistr->Add(hcosthpointd0d0S);
448     fDistr->Add(hcosthpointd0d0B);
449
450
451     //histograms of invariant mass distributions
452
453     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
454     TH1F *tmpMl=(TH1F*)tmpMt->Clone();
455     tmpMt->Sumw2();
456     tmpMl->Sumw2();
457
458     //to compare with AliAnalysisTaskCharmFraction
459     TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",200,1.765,1.965);
460     TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
461     tmpS27t->Sumw2();
462     tmpS27l->Sumw2();
463  
464     //distribution w/o cuts
465     //   TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
466     TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.56484,2.16484); //range (MD0-300MeV, mD0 + 300MeV)
467     TH1F *tmpMB=(TH1F*)tmpMS->Clone();
468     tmpMB->SetName(nameMassNocutsB.Data());
469     tmpMS->Sumw2();
470     tmpMB->Sumw2();
471
472     //MC signal and background
473     TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
474     TH1F *tmpSl=(TH1F*)tmpSt->Clone();
475     tmpSt->Sumw2();
476     tmpSl->Sumw2();
477
478     TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
479     TH1F *tmpBl=(TH1F*)tmpBt->Clone();
480     tmpBt->Sumw2();
481     tmpBl->Sumw2();
482
483     //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
484     TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
485     TH1F *tmpRl=(TH1F*)tmpRt->Clone();
486     tmpRt->Sumw2();
487     tmpRl->Sumw2();
488   //  printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
489
490     fOutputMass->Add(tmpMt);
491     fOutputMass->Add(tmpSt);
492     fOutputMass->Add(tmpS27t);
493     fOutputMass->Add(tmpBt);
494     fOutputMass->Add(tmpRt);
495
496     fDistr->Add(tmpMS);
497     fDistr->Add(tmpMB);
498
499
500   }
501
502   namedistr="hpospair";
503   TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
504   namedistr="hnegpair";
505   TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
506   fDistr->Add(hpospair);
507   fDistr->Add(hnegpair);
508
509
510   //histograms for vertex checking and TOF checking
511   TString checkname="hptGoodTr";
512
513   TH1F* hptGoodTr=new TH1F(checkname.Data(),"Pt distribution of 'good' tracks;p_{t}[GeV];Number",200,0.,8.);
514   hptGoodTr->SetTitleOffset(1.3,"Y");
515   checkname="hdistrGoodTr";
516
517   TH1F* hdistrGoodTr=new TH1F(checkname.Data(),"Distribution of number of good tracks per event;no.good-tracks/ev;Entries",31,0,31);
518   hdistrGoodTr->SetTitleOffset(1.3,"Y");
519
520   checkname="hTOFsig";
521   TH1F* hTOFsig=new TH1F(checkname.Data(),"Distribution of TOF signal;TOF time [ps];Entries", 100, -2.e3,40.e3);
522
523   checkname="hTPCsig";
524   TH1F* hTPCsig=new TH1F(checkname.Data(),"Distribution of TPC signal;TPC sig;Entries", 100, 35.,100.);
525
526   checkname="hTOFtimeKaonHyptime";
527   TH2F* hTOFtimeKaonHyptime=new TH2F(checkname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
528
529   checkname="hTOFtimeKaonHyptimeAC";
530   TH2F* hTOFtimeKaonHyptimeAC=new TH2F(checkname.Data(),"TOFtime - timeHypothesisForKaon;p[GeV/c];TOFtime - timeHypothesisForKaon [ps]",200,0.,4.,1000,-20000.,20000.);
531
532   checkname="hTPCsigvsp";
533   TH2F* hTPCsigvsp=new TH2F(checkname.Data(),"TPCsig vs p;p[GeV/c];TPCsig",200,0.,4.,1000,35.,100.);
534  
535  checkname="hTOFtime";
536   TH1F* hTOFtime=new TH1F(checkname.Data(),"Distribution of TOF time Kaon;TOF time(Kaon) [ps];Entries", 1000, 0.,50000.);
537   
538
539   fChecks->Add(hptGoodTr);
540   fChecks->Add(hdistrGoodTr);
541   fChecks->Add(hTOFsig);
542   fChecks->Add(hTPCsig);
543   fChecks->Add(hTOFtimeKaonHyptime);
544   fChecks->Add(hTOFtimeKaonHyptimeAC);
545   fChecks->Add(hTPCsigvsp);
546   fChecks->Add(hTOFtime);
547
548   const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
549
550   fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex ***  Integral(5,6) = pt out of bounds", 14,-0.5,13.5);
551
552   fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
553   fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
554   fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
555   fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
556   fNentries->GetXaxis()->SetBinLabel(5,"nEventsGoodVtx+>2tracks");
557   fNentries->GetXaxis()->SetBinLabel(6,"ptbin = -1");
558   fNentries->GetXaxis()->SetBinLabel(7,"no daughter");
559   fNentries->GetXaxis()->SetBinLabel(8,"nCandSel(Tr)");
560   fNentries->GetXaxis()->SetBinLabel(9,"PID=0");
561   fNentries->GetXaxis()->SetBinLabel(10,"PID=1");
562   fNentries->GetXaxis()->SetBinLabel(11,"PID=2");
563   fNentries->GetXaxis()->SetBinLabel(12,"PID=3");
564   fNentries->GetXaxis()->SetBinLabel(13,"K");
565   fNentries->GetXaxis()->SetBinLabel(14,"Lambda");
566   fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
567
568   fCounter = new AliNormalizationCounter("NormalizationCounter");
569
570   // Post the data
571   PostData(1,fOutputMass);
572   PostData(2,fDistr);
573   PostData(3,fNentries);
574   PostData(4,fChecks);
575   PostData(6,fCounter);  
576   return;
577 }
578
579 //________________________________________________________________________
580 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
581 {
582   // Execute analysis for current event:
583   // heavy flavor candidates association to MC truth
584   //cout<<"I'm in UserExec"<<endl;
585
586
587     //cuts order
588     //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
589     //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
590     //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
591     //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
592     //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
593     //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
594     //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
595     //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
596     //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
597   
598
599   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
600
601   TString bname;
602   if(fArray==0){ //D0 candidates
603     // load D0->Kpi candidates
604     //cout<<"D0 candidates"<<endl;
605     bname="D0toKpi";
606   } else { //LikeSign candidates
607     // load like sign candidates
608     //cout<<"LS candidates"<<endl;
609     bname="LikeSign2Prong";
610   }
611
612   TClonesArray *inputArray=0;
613  
614   if(!aod && AODEvent() && IsStandardAOD()) {
615     // In case there is an AOD handler writing a standard AOD, use the AOD 
616     // event in memory rather than the input (ESD) event.    
617     aod = dynamic_cast<AliAODEvent*> (AODEvent());
618     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
619     // have to taken from the AOD event hold by the AliAODExtension
620     AliAODHandler* aodHandler = (AliAODHandler*) 
621       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
622
623     if(aodHandler->GetExtensions()) {
624       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
625       AliAODEvent* aodFromExt = ext->GetAOD();
626       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
627     }
628   } else {
629     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
630   }
631
632
633   if(!inputArray) {
634     printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
635     return;
636   }
637   
638   // fix for temporary bug in ESDfilter
639   // the AODs with null vertex pointer didn't pass the PhysSel
640   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
641
642
643   TClonesArray *mcArray = 0;
644   AliAODMCHeader *mcHeader = 0;
645
646   if(fReadMC) {
647     // load MC particles
648     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
649     if(!mcArray) {
650       printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
651       return;
652     }
653     
654     // load MC header
655     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
656     if(!mcHeader) {
657       printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
658       return;
659     }
660   }
661   
662   //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
663   
664   //histogram filled with 1 for every AOD
665   fNentries->Fill(0);
666   fCounter->StoreEvent(aod,fReadMC);   
667   if(!fCuts->IsEventSelected(aod)) return;
668   
669   // AOD primary vertex
670   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
671
672   Bool_t isGoodVtx=kFALSE;
673
674    //vtx1->Print();
675   TString primTitle = vtx1->GetTitle();
676   if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
677     isGoodVtx=kTRUE;
678     fNentries->Fill(3);
679   }
680
681   
682   //cout<<"Start checks"<<endl;
683   Int_t ntracks=0,isGoodTrack=0;
684
685   if(aod) ntracks=aod->GetNTracks();
686
687   //cout<<"ntracks = "<<ntracks<<endl;
688  
689   //loop on tracks in the event
690   for (Int_t k=0;k<ntracks;k++){
691     AliAODTrack* track=aod->GetTrack(k);
692     //cout<<"track n = "<<k<<endl;
693  
694     //check TOF
695     if(!(track->GetStatus()&AliESDtrack::kTPCrefit &&
696          track->GetStatus()&AliESDtrack::kITSrefit && 
697          track->GetTPCNcls() >=70 &&
698          track->GetStatus()&AliESDtrack::kTOFpid && 
699          track->GetStatus()&AliESDtrack::kTOFout &&
700          track->GetStatus()&AliESDtrack::kTIME)) continue;
701     AliAODPid *pid = track->GetDetPid();
702     if(!pid)  {if (fDebug>1)cout<<"No AliAODPid found"<<endl; continue;}
703
704     Double_t times[5];
705     pid->GetIntegratedTimes(times);
706
707     ((TH1F*)fChecks->FindObject("hTOFtime"))->Fill(times[3]);
708     ((TH1F*)fChecks->FindObject("hTOFtimeKaonHyptime"))->Fill(track->P(),pid->GetTOFsignal()-times[3]); //3 is kaon
709
710     ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(pid->GetTOFsignal());
711     ((TH1F*)fChecks->FindObject("hTPCsig"))->Fill(pid->GetTPCsignal());
712     ((TH1F*)fChecks->FindObject("hTPCsigvsp"))->Fill(track->P(),pid->GetTPCsignal());
713     if (pid->GetTOFsignal()< 0) ((TH1F*)fChecks->FindObject("hTOFsig"))->Fill(-1);
714
715     //check clusters of the tracks
716     Int_t nclsTot=0,nclsSPD=0;
717     
718     for(Int_t l=0;l<6;l++) {
719       if(TESTBIT(track->GetITSClusterMap(),l)) {
720         nclsTot++; if(l<2) nclsSPD++;
721       }
722     }
723
724     if (track->Pt()>0.3 &&
725         track->GetStatus()&AliESDtrack::kTPCrefit &&
726         track->GetStatus()&AliESDtrack::kITSrefit &&
727         nclsTot>3 &&
728         nclsSPD>0) {//fill hist good tracks
729      
730       ((TH1F*)fChecks->FindObject("hptGoodTr"))->Fill(track->Pt());
731       isGoodTrack++;
732     }
733     //cout<<"isGoodTrack = "<<isGoodTrack<<endl;
734     ((TH1F*)fChecks->FindObject("hdistrGoodTr"))->Fill(isGoodTrack);
735   }
736   //number of events with good vertex and at least 2 good tracks
737   if (isGoodTrack>=2 && isGoodVtx) fNentries->Fill(4);
738
739   // loop over candidates
740   Int_t nInD0toKpi = inputArray->GetEntriesFast();
741   if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
742
743   // FILE *f=fopen("4display.txt","a");
744   // fprintf(f,"Number of D0->Kpi: %d\n",nInD0toKpi);
745   Int_t nSelectedloose=0,nSelectedtight=0;  
746   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
747     //Int_t nPosPairs=0, nNegPairs=0;
748     //cout<<"inside the loop"<<endl;
749     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
750
751     //check daughters
752     if(!(d->GetDaughter(0) || d->GetDaughter(1))) {
753       AliDebug(1,"at least one daughter not found!");
754       fNentries->Fill(6);
755       continue;
756     }
757
758     // Bool_t unsetvtx=kFALSE;
759     // if(!d->GetOwnPrimaryVtx()) {
760     //   d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
761     //   unsetvtx=kTRUE;
762     // }
763   
764     
765     //check reco daughter in acceptance
766     /*
767     Double_t eta0=d->EtaProng(0);
768     Double_t eta1=d->EtaProng(1);
769     */
770     if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
771       //if( TMath::Abs(eta0)<0.9 && TMath::Abs(eta1)<0.9 ){
772        //apply cuts on tracks
773       nSelectedloose++;
774       Int_t isSelected = fCuts->IsSelected(d,AliRDHFCuts::kTracks);
775       if(((AliAODTrack*)d->GetDaughter(0))->GetTPCNcls() < 70 || ((AliAODTrack*)d->GetDaughter(1))->GetTPCNcls() < 70) isSelected=kFALSE;
776       if (!isSelected) continue;
777       nSelectedtight++;
778       fNentries->Fill(7);       
779       if(fDebug>2) cout<<"tracks selected"<<endl;
780
781       Int_t ptbin=fCuts->PtBin(d->Pt());
782       if(ptbin==-1) {fNentries->Fill(5); continue;} //out of bounds
783       FillVarHists(aod,d,mcArray,fCuts,fDistr);
784       FillMassHists(aod,d,mcArray,fCuts,fOutputMass);
785     }
786   
787     //if(unsetvtx) d->UnsetOwnPrimaryVtx();
788   } //end for prongs
789   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);  
790   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);  
791   // Post the data
792   PostData(1,fOutputMass);
793   PostData(2,fDistr);
794   PostData(3,fNentries);
795   PostData(4,fChecks);
796   PostData(6,fCounter);
797   return;
798 }
799
800 //____________________________________________________________________________
801 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
802   //
803   // function used in UserExec to fill variable histograms:
804   //
805
806   //add distr here
807   UInt_t pdgs[2];
808     
809   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
810   pdgs[0]=211;
811   pdgs[1]=321;
812   Double_t minvD0 = part->InvMassD0();
813   pdgs[1]=211;
814   pdgs[0]=321;
815   Double_t minvD0bar = part->InvMassD0bar();
816   //cout<<"inside mass cut"<<endl;
817
818   Int_t pdgDgD0toKpi[2]={321,211};
819   Int_t lab=-9999;
820   if(fReadMC) lab=part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
821   //Double_t pt = d->Pt(); //mother pt
822   Int_t isSelected=3;
823   Int_t isSelectedPID=99;
824
825   //isSelectedPID = cuts->IsSelected(part,AliRDHFCuts::kPID); //0 rejected,1 D0,2 Dobar, 3 both
826   isSelectedPID = cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
827   if (isSelectedPID==0)fNentries->Fill(8);
828   if (isSelectedPID==1)fNentries->Fill(9);
829   if (isSelectedPID==2)fNentries->Fill(10);
830   if (isSelectedPID==3)fNentries->Fill(11);
831     //fNentries->Fill(8+isSelectedPID);
832
833   if(fCutOnDistr){
834     isSelected = cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod); //cuts with variables recalculated with new vertex (w/o daughters)
835     if (!isSelected){
836       //cout<<"Not Selected"<<endl;
837       return;
838     }
839   }
840
841   Double_t invmasscut=0.03;
842
843   TString fillthispi="",fillthisK="",fillthis="";
844
845   Int_t ptbin=cuts->PtBin(part->Pt());
846
847
848   //recalculate vertex
849   AliAODVertex *vtxProp=0x0;
850   vtxProp=GetPrimaryVtxSkipped(aod,part);
851   Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
852   dz1[0]=-99; dz2[0]=-99;
853   if(vtxProp) {
854     part->SetOwnPrimaryVtx(vtxProp);
855     //Bool_t unsetvtx=kTRUE;
856     //Calculate d0 for daughter tracks with Vtx Skipped
857     AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)part->GetDaughter(0));
858     AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)part->GetDaughter(1));
859     esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
860     esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
861     delete vtxProp; vtxProp=NULL;
862     delete esdtrack1;
863     delete esdtrack2;
864   }
865
866   Double_t d0[2]={dz1[0],dz2[0]};
867   Double_t decl[2]={part->DecayLength(),part->NormalizedDecayLength()};
868   part->UnsetOwnPrimaryVtx();
869
870   if(!fCutOnDistr || (fCutOnDistr && isSelected)){ //if no cuts or cuts passed 
871     //printf("\nif no cuts or cuts passed\n");
872
873     //disable the PID
874     if(!fUsePid4Distr) isSelectedPID=0;
875     if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
876       //if(!vtxProp) fNentries->Fill(3);
877       //check pdg of the prongs
878       AliAODTrack *prong0=(AliAODTrack*)part->GetDaughter(0);
879       AliAODTrack *prong1=(AliAODTrack*)part->GetDaughter(1);
880       if(!prong0 || !prong1) {
881         return;
882       }
883       Int_t labprong[2];
884       if(fReadMC){
885         labprong[0]=prong0->GetLabel();
886         labprong[1]=prong1->GetLabel();
887       }
888       AliAODMCParticle *mcprong=0;
889       Int_t pdgProng[2]={0,0};
890
891       for (Int_t iprong=0;iprong<2;iprong++){
892         if(fReadMC && labprong[iprong]>=0) {
893           mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
894           pdgProng[iprong]=mcprong->GetPdgCode();
895         }
896       }
897     
898
899       //no mass cut ditributions: ptbis
900         
901       fillthispi="hptpiSnoMcut_";
902       fillthispi+=ptbin;
903
904       fillthisK="hptKSnoMcut_";
905       fillthisK+=ptbin;
906
907       if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
908           || (isSelectedPID==1 || isSelectedPID==3)){
909         ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(0));
910         ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(1));
911       }
912
913       if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
914           || (isSelectedPID==2 || isSelectedPID==3)){
915         ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(0));
916         ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(1));
917       }
918       
919       //no mass cut ditributions: mass
920       fillthis="hMassS_";
921       fillthis+=ptbin;
922       
923       if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
924           || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
925         ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
926       }
927       else { //D0bar
928         if(fReadMC || (!fReadMC && isSelectedPID > 1))
929         ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
930       }
931
932       //apply cut on invariant mass on the pair
933       if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
934         
935         if(fArray==1) cout<<"LS signal: ERROR"<<endl;
936         for (Int_t iprong=0; iprong<2; iprong++){
937           AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
938           if (fReadMC) labprong[iprong]=prong->GetLabel();
939           
940           //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
941           Int_t pdgprong=0;
942           if(fReadMC && labprong[iprong]>=0) {
943             mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
944             pdgprong=mcprong->GetPdgCode();
945           }
946
947           Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
948
949           if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
950             //cout<<"pi"<<endl;
951             
952             fillthispi="hptpiS_";
953             fillthispi+=ptbin;
954             ((TH1F*)listout->FindObject(fillthispi))->Fill(part->PtProng(iprong));
955             fillthispi="hd0piS_";
956             fillthispi+=ptbin;
957             ((TH1F*)listout->FindObject(fillthispi))->Fill(part->Getd0Prong(iprong));
958             if(d0[iprong] != -99) {
959
960               fillthispi="hd0vpiS_";
961               fillthispi+=ptbin;
962               ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
963             }
964
965           }
966           
967           if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
968             //cout<<"kappa"<<endl;
969             
970             fillthisK="hptKS_";
971             fillthisK+=ptbin;
972             ((TH1F*)listout->FindObject(fillthisK))->Fill(part->PtProng(iprong));
973             fillthisK="hd0KS_";
974             fillthisK+=ptbin;
975             ((TH1F*)listout->FindObject(fillthisK))->Fill(part->Getd0Prong(iprong));
976             if (d0[iprong] != -99){
977               fillthisK="hd0vKS_";
978               fillthisK+=ptbin;
979               ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
980             }
981           }
982
983           fillthis="hcosthpointd0S_";
984           fillthis+=ptbin;        
985           ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(iprong));
986
987         } //end loop on prongs
988
989         fillthis="hdcaS_";
990         fillthis+=ptbin;          
991         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
992
993         fillthis="hcosthetapointS_";
994         fillthis+=ptbin;          
995         ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
996
997         fillthis="hcosthpointd0d0S_";
998         fillthis+=ptbin;          
999         ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
1000           
1001         fillthis="hcosthetastarS_";
1002         fillthis+=ptbin;
1003         if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
1004         else {
1005           if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());
1006           if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
1007         }
1008         fillthis="hd0d0S_";
1009         fillthis+=ptbin;
1010         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1011
1012         if(d0[0] != -99 && d0[1] != -99 ){
1013           fillthis="hd0d0vS_";
1014           fillthis+=ptbin;
1015           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1016         }
1017
1018         fillthis="hdeclS_";
1019         fillthis+=ptbin;
1020         ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
1021
1022         fillthis="hnormdeclS_";
1023         fillthis+=ptbin;
1024         ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
1025
1026         if(vtxProp) {
1027           fillthis="hdeclvS_";
1028           fillthis+=ptbin;
1029           ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1030
1031           fillthis="hnormdeclvS_";
1032           fillthis+=ptbin;
1033           ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1034         }
1035       } //end mass cut
1036     
1037     } else{ //Background or LS
1038       //if(!fReadMC){
1039       //cout<<"is background"<<endl;
1040      
1041       //no mass cut distributions: mass, ptbis
1042       fillthis="hMassB_";
1043       fillthis+=ptbin;
1044       
1045       if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3))) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1046       if (!fCutOnDistr || (fCutOnDistr && isSelected>1)) ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1047
1048       fillthis="hptB1prongnoMcut_";
1049       fillthis+=ptbin;
1050       
1051       ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
1052       
1053       fillthis="hptB2prongsnoMcut_";
1054       fillthis+=ptbin;
1055       ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
1056       ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
1057
1058       //apply cut on invariant mass on the pair
1059       if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1060
1061
1062         AliAODTrack *prongg=(AliAODTrack*)part->GetDaughter(0);
1063         if(!prongg) {
1064           if(fDebug>2) cout<<"No daughter found";
1065           return;
1066         }
1067         else{
1068           if(prongg->Charge()==1) {
1069             fTotPosPairs[ptbin]++;
1070             ((TH1F*)fDistr->FindObject("hpospair"))->Fill(ptbin);
1071           } else {
1072             fTotNegPairs[ptbin]++;
1073             ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(ptbin);
1074           }
1075         }
1076         
1077         //normalise pt distr to half afterwards
1078         fillthis="hptB_";
1079         fillthis+=ptbin;
1080         ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(0));
1081         ((TH1F*)listout->FindObject(fillthis))->Fill(part->PtProng(1));
1082
1083         fillthis="hd0p0B_";
1084         fillthis+=ptbin;
1085         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
1086         fillthis="hd0p1B_";
1087         fillthis+=ptbin;
1088         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
1089         fillthis="hd0B_";
1090         fillthis+=ptbin;
1091         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(0));
1092         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(1));
1093
1094         if(fReadMC){
1095           Int_t pdgMother[2]={0,0};
1096           Double_t factor[2]={1,1};
1097
1098           for(Int_t iprong=0;iprong<2;iprong++){
1099             AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(iprong);
1100             lab=prong->GetLabel();
1101             if(lab>=0){
1102               AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
1103               if(mcprong){
1104                 Int_t labmom=mcprong->GetMother();
1105                 if(labmom>=0){
1106                   AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
1107                   if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
1108                 }
1109               }
1110             }
1111
1112             fillthis="hd0moresB_";
1113             fillthis+=ptbin;
1114           
1115             if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
1116               if(part->PtProng(iprong)<=1)factor[iprong]=1./.7;
1117               else factor[iprong]=1./.6;
1118               fNentries->Fill(12);
1119             }
1120             
1121             if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
1122               factor[iprong]=1./0.25;
1123               fNentries->Fill(13);
1124             }
1125             fillthis="hd0moresB_";
1126             fillthis+=ptbin;
1127
1128             ((TH1F*)listout->FindObject(fillthis))->Fill(part->Getd0Prong(iprong),factor[iprong]);
1129
1130             fillthis="hd0vmoresB_";
1131             fillthis+=ptbin;
1132             ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
1133
1134           }
1135
1136           fillthis="hd0d0moresB_";
1137           fillthis+=ptbin;
1138           ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1139
1140           fillthis="hd0d0vmoresB_";
1141           fillthis+=ptbin;
1142           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1143         
1144           fillthis="hcosthetapointmoresB_";
1145           fillthis+=ptbin;
1146           ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),factor[0]*factor[1]);
1147         }
1148
1149         fillthis="hd0vp0B_";
1150         fillthis+=ptbin;
1151         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1152         fillthis="hd0vp1B_";
1153         fillthis+=ptbin;
1154         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1155
1156         fillthis="hd0vB_";
1157         fillthis+=ptbin;
1158         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1159         ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1160   
1161
1162         fillthis="hcosthpointd0B_";
1163         fillthis+=ptbin;          
1164         ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(0));
1165         ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Getd0Prong(1));
1166
1167
1168         fillthis="hdcaB_";
1169         fillthis+=ptbin;
1170         ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1171
1172         fillthis="hcosthetastarB_";
1173         fillthis+=ptbin;
1174         if (!fCutOnDistr || (fCutOnDistr && (isSelected==1 || isSelected==3)))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0());
1175         if (!fCutOnDistr || (fCutOnDistr && isSelected>1))((TH1F*)listout->FindObject(fillthis))->Fill(part->CosThetaStarD0bar());      
1176
1177         fillthis="hd0d0B_";
1178         fillthis+=ptbin;
1179         ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1180
1181         if(d0[0] != -99 && d0[1] != -99 ){
1182           fillthis="hd0d0vB_";
1183           fillthis+=ptbin;
1184           ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1185         }
1186
1187         fillthis="hcosthetapointB_";
1188         fillthis+=ptbin;
1189         ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle());
1190
1191         fillthis="hcosthpointd0d0B_";
1192         fillthis+=ptbin;
1193         ((TH2F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngle(),part->Prodd0d0());
1194
1195         fillthis="hdeclB_";
1196         fillthis+=ptbin;
1197         ((TH1F*)listout->FindObject(fillthis))->Fill(part->DecayLength());
1198
1199         fillthis="hnormdeclB_";
1200         fillthis+=ptbin;
1201         ((TH1F*)listout->FindObject(fillthis))->Fill(part->NormalizedDecayLength());
1202
1203         if(vtxProp) {
1204           fillthis="hdeclvB_";
1205           fillthis+=ptbin;
1206           ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1207
1208           fillthis="hnormdeclvB_";
1209           fillthis+=ptbin;
1210           ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1211         }
1212       }//mass cut       
1213     }//else (background)
1214   }
1215   return;
1216 }
1217 //____________________________________________________________________________
1218
1219 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi* cuts, TList *listout){
1220   //
1221   // function used in UserExec to fill mass histograms:
1222   //
1223
1224
1225   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1226
1227   Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate,aod); //selected
1228
1229   //cout<<"is selected = "<<isSelected<<endl;
1230
1231   //cout<<"check cuts = "<<endl;
1232   //cuts->PrintAll();
1233   if (!isSelected){
1234     //cout<<"Not Selected"<<endl;
1235     //cout<<"Rejected because "<<cuts->GetWhy()<<endl;
1236     return;
1237   }
1238
1239
1240   if(fDebug>2)  cout<<"Candidate selected"<<endl;
1241
1242   Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
1243   //printf("SELECTED\n");
1244   Int_t ptbin=cuts->PtBin(part->Pt());
1245
1246   AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
1247   if(!prong) {
1248     AliDebug(2,"No daughter found");
1249     return;
1250   }
1251   else{
1252     if(prong->Charge()==1) {fTotPosPairs[ptbin]++;} else {fTotNegPairs[ptbin]++;}
1253   }
1254
1255   for(Int_t it=0;it<2;it++){
1256     //fill TOFsig-kaontime
1257     AliAODTrack* track=(AliAODTrack*)part->GetDaughter(it);
1258     AliAODPid *pid = track->GetDetPid();
1259     Double_t times[5];
1260     pid->GetIntegratedTimes(times);
1261     ((TH2F*)fChecks->FindObject("hTOFtimeKaonHyptimeAC"))->Fill(track->P(),pid->GetTOFsignal()-times[3]);
1262   }
1263
1264   //request on spd points to be addes
1265   if(/*nSPD==2 && */part->Pt() > 5. && (TMath::Abs(invmassD0-mPDG)<0.01 || TMath::Abs(invmassD0bar-mPDG)<0.01)){
1266     FILE *f=fopen("4display.txt","a");
1267     fprintf(f,"pt: %f \n Rapidity: %f \t Period Number: %x \t Run Number: %d \t BunchCrossNumb: %d \t OrbitNumber: %d\n",part->Pt(),part->Y(421),aod->GetPeriodNumber(),aod->GetRunNumber(),aod->GetBunchCrossNumber(),aod->GetOrbitNumber());
1268     fclose(f);
1269     //printf("PrimVtx NContributors: %d \n Prongs Rel Angle: %f \n \n",ncont,relangle);
1270   }
1271  
1272
1273
1274   TString fillthis="";
1275   Int_t pdgDgD0toKpi[2]={321,211};
1276   Int_t labD0=-1;
1277   if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
1278
1279   //count candidates selected by cuts
1280   fNentries->Fill(1);
1281   //count true D0 selected by cuts
1282   if (fReadMC && labD0>=0) fNentries->Fill(2);
1283   //PostData(3,fNentries);
1284
1285   if ((isSelected==1 || isSelected==3) && fFillOnlyD0D0bar<2) { //D0
1286     fillthis="histMass_";
1287     fillthis+=ptbin;
1288     //cout<<"Filling "<<fillthis<<endl;
1289
1290     //printf("Fill mass with D0");
1291     ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1292  
1293     if(labD0>=0) {
1294       if(fArray==1) cout<<"LS signal ERROR"<<endl;
1295
1296       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1297       Int_t pdgD0 = partD0->GetPdgCode();
1298       //cout<<"pdg = "<<pdgD0<<endl;
1299       if (pdgD0==421){ //D0
1300         //cout<<"Fill S with D0"<<endl;
1301         fillthis="histSgn_";
1302         fillthis+=ptbin;
1303         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1304         if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1305           fillthis="histSgn27_";
1306           fillthis+=ptbin;
1307           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1308         }
1309       } else{ //it was a D0bar
1310         fillthis="histRfl_";
1311         fillthis+=ptbin;
1312         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1313       }
1314     } else {//background
1315       fillthis="histBkg_";
1316       fillthis+=ptbin;
1317       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1318     }
1319       
1320   }
1321   if (isSelected>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
1322     fillthis="histMass_";
1323     fillthis+=ptbin;
1324     //printf("Fill mass with D0bar");
1325     ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
1326  
1327     if(labD0>=0) {
1328       if(fArray==1) cout<<"LS signal ERROR"<<endl;
1329       AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1330       Int_t pdgD0 = partD0->GetPdgCode();
1331       //cout<<" pdg = "<<pdgD0<<endl;
1332       if (pdgD0==-421){ //D0bar
1333         fillthis="histSgn_";
1334         fillthis+=ptbin;
1335         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1336         if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1337           fillthis="histSgn27_";
1338           fillthis+=ptbin;
1339           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1340         }
1341
1342           
1343       } else{
1344         fillthis="histRfl_";
1345         fillthis+=ptbin;
1346         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1347       }
1348     } else {//background or LS
1349       fillthis="histBkg_";
1350       fillthis+=ptbin;
1351       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1352     }
1353   }
1354
1355   return;
1356 }
1357
1358 //__________________________________________________________________________
1359 AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev,AliAODRecoDecayHF2Prong *d){
1360   //Calculate the primary vertex w/o the daughter tracks of the candidate
1361   
1362   Int_t skipped[2];
1363   Int_t nTrksToSkip=2;
1364   AliAODTrack *dgTrack = (AliAODTrack*)d->GetDaughter(0);
1365   if(!dgTrack){
1366     AliDebug(2,"no daughter found!");
1367     return 0x0;
1368   }
1369   skipped[0]=dgTrack->GetID();
1370   dgTrack = (AliAODTrack*)d->GetDaughter(1);
1371   if(!dgTrack){
1372     AliDebug(2,"no daughter found!");
1373     return 0x0;
1374   }
1375   skipped[1]=dgTrack->GetID();
1376
1377   AliESDVertex *vertexESD=0x0;
1378   AliAODVertex *vertexAOD=0x0;
1379   AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
1380   
1381   //
1382   vertexer->SetSkipTracks(nTrksToSkip,skipped);
1383   vertexer->SetMinClusters(4);  
1384   vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev); 
1385   if(!vertexESD) return vertexAOD;
1386   if(vertexESD->GetNContributors()<=0) { 
1387     AliDebug(2,"vertexing failed"); 
1388     delete vertexESD; vertexESD=NULL;
1389     return vertexAOD;
1390   }
1391   
1392   delete vertexer; vertexer=NULL;
1393   
1394   
1395   // convert to AliAODVertex
1396   Double_t pos[3],cov[6],chi2perNDF;
1397   vertexESD->GetXYZ(pos); // position
1398   vertexESD->GetCovMatrix(cov); //covariance matrix
1399   chi2perNDF = vertexESD->GetChi2toNDF();
1400   delete vertexESD; vertexESD=NULL;
1401   
1402   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
1403   return vertexAOD;
1404   
1405 }
1406
1407
1408 //________________________________________________________________________
1409 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1410 {
1411   // Terminate analysis
1412   //
1413   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1414
1415
1416   fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
1417   if (!fOutputMass) {     
1418     printf("ERROR: fOutputMass not available\n");
1419     return;
1420   }
1421   fDistr = dynamic_cast<TList*> (GetOutputData(2));
1422   if (!fDistr) {
1423     printf("ERROR: fDistr not available\n");
1424     return;
1425   }
1426  
1427   fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
1428   
1429   if(!fNentries){
1430     printf("ERROR: fNEntries not available\n");
1431     return;
1432   }
1433
1434   fChecks = dynamic_cast<TList*> (GetOutputData(4));
1435   if (!fChecks) {
1436     printf("ERROR: fChecks not available\n");
1437     return;
1438   }
1439   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(6));    
1440   if (!fCounter) {
1441     printf("ERROR: fCounter not available\n");
1442     return;
1443   }
1444   for(Int_t ipt=0;ipt<5;ipt++){ //change 5 in GetNPtBins when sure it is written and check
1445
1446
1447     if(fArray==1){ 
1448       fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]); 
1449
1450
1451       if(fLsNormalization>1e-6) {
1452         
1453         TString massName="histMass_";
1454         massName+=ipt;
1455         ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
1456
1457       }
1458     
1459
1460       fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1461
1462       if(fLsNormalization>1e-6) {
1463
1464         TString nameDistr="hptB_";
1465         nameDistr+=ipt;
1466         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1467         nameDistr="hdcaB_";
1468         nameDistr+=ipt;
1469         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1470         nameDistr="hcosthetastarB_";
1471         nameDistr+=ipt;
1472         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1473         nameDistr="hd0B_";
1474         nameDistr+=ipt;
1475         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1476         nameDistr="hd0d0B_";
1477         nameDistr+=ipt;
1478         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1479         nameDistr="hcosthetapointB_";
1480         nameDistr+=ipt;
1481         ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1482         nameDistr="hcosthpointd0d0B_";
1483         nameDistr+=ipt;
1484         ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
1485
1486       }
1487     }
1488   }
1489   TString cvname,cstname;
1490
1491   if (fArray==0){
1492     cvname="D0invmass";
1493     cstname="cstat0";
1494   } else {
1495     cvname="LSinvmass";
1496     cstname="cstat1";
1497   }
1498
1499   TCanvas *cMass=new TCanvas(cvname,cvname);
1500   cMass->cd();
1501   ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
1502
1503   TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
1504   cStat->cd();
1505   cStat->SetGridy();
1506   fNentries->Draw("htext0");
1507
1508   // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
1509   // ccheck->cd();
1510
1511   return;
1512 }
1513