Updates (Chiara)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSED0Mass.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /////////////////////////////////////////////////////////////
18 //
19 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
20 // and comparison with the MC truth and cut variables distributions.
21 //
22 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
23 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
24 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
25 /////////////////////////////////////////////////////////////
26
27 #include <Riostream.h>
28 #include <TClonesArray.h>
29 #include <TCanvas.h>
30 #include <TNtuple.h>
31 #include <TList.h>
32 #include <TH1F.h>
33 #include <TH2F.h>
34 #include <TDatabasePDG.h>
35
36 #include "AliAnalysisManager.h"
37 #include "AliAODHandler.h"
38 #include "AliAODEvent.h"
39 #include "AliAODVertex.h"
40 #include "AliAODTrack.h"
41 #include "AliAODMCHeader.h"
42 #include "AliAODMCParticle.h"
43 #include "AliAODRecoDecayHF2Prong.h"
44 #include "AliAODRecoCascadeHF.h"
45 #include "AliAnalysisVertexingHF.h"
46 #include "AliAnalysisTaskSE.h"
47 #include "AliAnalysisTaskSED0Mass.h"
48
49
50 ClassImp(AliAnalysisTaskSED0Mass)
51
52
53 //________________________________________________________________________
54 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
55 AliAnalysisTaskSE(),
56 fOutputPPR(0),
57 fOutputmycuts(0),
58 fDistr(0), 
59 fNentries(0),
60 fVHFPPR(0),
61 fVHFmycuts(0),
62 fArray(0),
63 fReadMC(0),
64 fLsNormalization(1.)
65
66 {
67   // Default constructor
68   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
69 }
70
71 //________________________________________________________________________
72 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name):
73 AliAnalysisTaskSE(name),
74 fOutputPPR(0), 
75 fOutputmycuts(0), 
76 fDistr(0),
77 fNentries(0),
78 fVHFPPR(0),
79 fVHFmycuts(0),
80 fArray(0),
81 fReadMC(0),
82 fLsNormalization(1.)
83
84 {
85   // Default constructor
86   for(Int_t i=0;i<5;i++) {fTotPosPairs[i]=0; fTotNegPairs[i]=0;}
87
88   // Output slot #1 writes into a TList container
89   DefineOutput(1,TList::Class());  //My private output
90   // Output slot #2 writes into a TList container
91   DefineOutput(2,TList::Class());  //My private output
92   // Output slot #3 writes into a TH1F container
93   DefineOutput(3,TH1F::Class());  //My private output
94   // Output slot #4 writes into a TList container
95   DefineOutput(4,TList::Class());  //My private output
96 }
97
98 //________________________________________________________________________
99 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
100 {
101    if (fOutputPPR) {
102     delete fOutputPPR;
103     fOutputPPR = 0;
104   }
105   if (fVHFPPR) {
106     delete fVHFPPR;
107     fVHFPPR = 0;
108   }
109   if (fOutputmycuts) {
110     delete fOutputmycuts;
111     fOutputmycuts = 0;
112   }
113
114   if (fDistr) {
115     delete fDistr;
116     fDistr = 0;
117   }
118
119   if (fVHFmycuts) {
120     delete fVHFmycuts;
121     fVHFmycuts = 0;
122   }
123   if (fNentries){
124     delete fNentries;
125     fNentries = 0;
126   }
127
128
129 }  
130
131 //________________________________________________________________________
132 void AliAnalysisTaskSED0Mass::Init()
133 {
134   // Initialization
135
136   if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
137
138   gROOT->LoadMacro("$ALICE_ROOT/PWG3/vertexingHF/ConfigVertexingHF.C");
139
140   // 2 sets of dedidcated cuts -- defined in UserExec
141   fVHFPPR = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
142   fVHFmycuts = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");  
143   
144   return;
145 }
146
147 //________________________________________________________________________
148 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
149 {
150
151   // Create the output container
152   //
153   if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
154
155   // Several histograms are more conveniently managed in a TList
156   fOutputPPR = new TList();
157   fOutputPPR->SetOwner();
158   fOutputPPR->SetName("listPPR");
159
160   fOutputmycuts = new TList();
161   fOutputmycuts->SetOwner();
162   fOutputmycuts->SetName("listloose");
163
164   fDistr = new TList();
165   fDistr->SetOwner();
166   fDistr->SetName("distributionslist");
167
168   const Int_t nhist=5;
169
170   TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
171
172   for(Int_t i=0;i<nhist;i++){
173
174     nameMass="histMass_";
175     nameMass+=i+1;
176     nameSgn27="histSgn27_";
177     nameSgn27+=i+1;
178     nameSgn="histSgn_";
179     nameSgn+=i+1;
180     nameBkg="histBkg_";
181     nameBkg+=i+1;
182     nameRfl="histRfl_";
183     nameRfl+=i+1;
184     nameMassNocutsS="hMassS_";
185     nameMassNocutsS+=i+1;
186     nameMassNocutsB="hMassB_";
187     nameMassNocutsB+=i+1;
188
189     //histograms of cut variable distributions
190
191     //  pT
192     namedistr="hptpiS_";
193     namedistr+=i+1;
194     TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
195
196     namedistr="hptKS_";
197     namedistr+=i+1;
198     TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
199
200     namedistr="hptB_";
201     namedistr+=i+1;
202     TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
203
204     //  pT no mass cut
205     namedistr="hptpiSnoMcut_";
206     namedistr+=i+1;
207     TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
208
209     namedistr="hptKSnoMcut_";
210     namedistr+=i+1;
211     TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
212
213     namedistr="hptB1prongnoMcut_";
214     namedistr+=i+1;
215     TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
216
217     namedistr="hptB2prongsnoMcut_";
218     namedistr+=i+1;
219     TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
220
221     
222     //  dca
223     namedistr="hdcaS_";
224     namedistr+=i+1;
225     TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
226     namedistr="hdcaB_";
227     namedistr+=i+1;
228     TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
229
230     //  costhetastar
231     namedistr="hcosthetastarS_";
232     namedistr+=i+1;
233     TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
234     namedistr="hcosthetastarB_";
235     namedistr+=i+1;
236     TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
237
238     // impact parameter
239     namedistr="hd0piS_";
240     namedistr+=i+1;
241     TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
242
243     namedistr="hd0KS_";
244     namedistr+=i+1;
245     TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
246     namedistr="hd0B_";
247     namedistr+=i+1;
248     TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution;d0 [cm]",200,-0.1,0.1);
249
250     namedistr="hd0d0S_";
251     namedistr+=i+1;
252     TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
253     namedistr="hd0d0B_";
254     namedistr+=i+1;
255     TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
256
257     //  costhetapoint
258     namedistr="hcosthetapointS_";
259     namedistr+=i+1;
260     TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
261     namedistr="hcosthetapointB_";
262     namedistr+=i+1;
263     TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
264
265     namedistr="hcosthpointd0d0S_";
266     namedistr+=i+1;
267     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);
268     namedistr="hcosthpointd0d0B_";
269     namedistr+=i+1;
270     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);
271
272     fDistr->Add(hptpiS);
273     fDistr->Add(hptKS);
274     fDistr->Add(hptB);
275
276     fDistr->Add(hptpiSnoMcut);
277     fDistr->Add(hptKSnoMcut);
278     fDistr->Add(hptB1pnoMcut);
279     fDistr->Add(hptB2pnoMcut);
280
281     fDistr->Add(hdcaS);
282     fDistr->Add(hdcaB);
283
284     fDistr->Add(hd0piS);
285     fDistr->Add(hd0KS);
286     fDistr->Add(hd0B);
287
288     fDistr->Add(hd0d0S);
289     fDistr->Add(hd0d0B);
290
291     fDistr->Add(hcosthetastarS);
292     fDistr->Add(hcosthetastarB);
293
294     fDistr->Add(hcosthetapointS);
295     fDistr->Add(hcosthetapointB);
296
297     fDistr->Add(hcosthpointd0d0S);
298     fDistr->Add(hcosthpointd0d0B);
299
300
301     //histograms of invariant mass distributions
302
303     TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",200,1.765,1.965);
304     TH1F *tmpMl=(TH1F*)tmpMt->Clone();
305     tmpMt->Sumw2();
306     tmpMl->Sumw2();
307
308     //to compare with AliAnalysisTaskCharmFraction
309     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);
310     TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
311     tmpS27t->Sumw2();
312     tmpS27l->Sumw2();
313  
314     //distribution w/o cuts
315     //   TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
316     TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.56484,2.16484); //range (MD0-300MeV, mD0 + 300MeV)
317     TH1F *tmpMB=(TH1F*)tmpMS->Clone();
318     tmpMB->SetName(nameMassNocutsB.Data());
319     tmpMS->Sumw2();
320     tmpMB->Sumw2();
321
322     //MC signal and background
323     TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
324     TH1F *tmpSl=(TH1F*)tmpSt->Clone();
325     tmpSt->Sumw2();
326     tmpSl->Sumw2();
327
328     TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
329     TH1F *tmpBl=(TH1F*)tmpBt->Clone();
330     tmpBt->Sumw2();
331     tmpBl->Sumw2();
332
333     //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
334     TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",200,1.765,1.965);
335     TH1F *tmpRl=(TH1F*)tmpRt->Clone();
336     tmpRt->Sumw2();
337     tmpRl->Sumw2();
338   //  printf("Created histograms %s\t%s\t%s\t%s\n",tmpM->GetName(),tmpS->GetName(),tmpB->GetName(),tmpR->GetName());
339
340     fOutputPPR->Add(tmpMt);
341     fOutputPPR->Add(tmpSt);
342     fOutputPPR->Add(tmpS27t);
343     fOutputPPR->Add(tmpBt);
344     fOutputPPR->Add(tmpRt);
345
346     fOutputmycuts->Add(tmpMl);
347     fOutputmycuts->Add(tmpSl);
348     fOutputmycuts->Add(tmpS27l);
349     fOutputmycuts->Add(tmpBl);
350     fOutputmycuts->Add(tmpRl);
351  
352     fDistr->Add(tmpMS);
353     fDistr->Add(tmpMB);
354
355   }
356
357
358   fNentries=new TH1F("nentriesD0", "nentriesD0->Integral(1,2) = number of AODs *** nentriesD0->Integral(3,4) = number of candidates selected with cuts *** nentriesD0->Integral(5,6) = number of D0 selected with cuts", 6,1.,4.);
359
360   return;
361 }
362
363 //________________________________________________________________________
364 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
365 {
366   // Execute analysis for current event:
367   // heavy flavor candidates association to MC truth
368   //cout<<"I'm in UserExec"<<endl;
369   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
370
371   TString bname;
372   if(fArray==0){ //D0 candidates
373     // load D0->Kpi candidates
374     //cout<<"D0 candidates"<<endl;
375     bname="D0toKpi";
376   } else { //LikeSign candidates
377     // load like sign candidates
378     //cout<<"LS candidates"<<endl;
379     bname="LikeSign2Prong";
380   }
381
382   TClonesArray *inputArray=0;
383
384   if(!aod && AODEvent() && IsStandardAOD()) {
385     // In case there is an AOD handler writing a standard AOD, use the AOD 
386     // event in memory rather than the input (ESD) event.    
387     aod = dynamic_cast<AliAODEvent*> (AODEvent());
388     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
389     // have to taken from the AOD event hold by the AliAODExtension
390     AliAODHandler* aodHandler = (AliAODHandler*) 
391       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
392     if(aodHandler->GetExtensions()) {
393       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
394       AliAODEvent *aodFromExt = ext->GetAOD();
395       inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
396     }
397   } else {
398     inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
399   }
400
401
402   if(!inputArray) {
403     printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
404     return;
405   }
406   
407   // AOD primary vertex
408   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
409   //vtx1->Print();
410
411   TClonesArray *mcArray = 0;
412   AliAODMCHeader *mcHeader = 0;
413
414   if(fReadMC) {
415     // load MC particles
416     mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
417     if(!mcArray) {
418       printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
419       return;
420     }
421   
422     // load MC header
423     mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
424     if(!mcHeader) {
425       printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
426       return;
427     }
428   }
429   
430   //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
431
432   //histogram filled with 1 for every AOD
433   fNentries->Fill(1);
434   PostData(3,fNentries);
435
436   // loop over candidates
437   Int_t nInD0toKpi = inputArray->GetEntriesFast();
438   if(fDebug>1) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
439   
440   for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
441     //Int_t nPosPairs=0, nNegPairs=0;
442     //cout<<"inside the loop"<<endl;
443     AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
444     Bool_t unsetvtx=kFALSE;
445     if(!d->GetOwnPrimaryVtx()) {
446       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
447       unsetvtx=kTRUE;
448     }
449     
450     //check reco daughter in acceptance
451     Double_t eta0=d->EtaProng(0);
452     Double_t eta1=d->EtaProng(1);
453     Bool_t prongsinacc=kFALSE;
454     if (TMath::Abs(eta0) < 0.9 && TMath::Abs(eta1) < 0.9) prongsinacc=kTRUE;
455
456     //add distr here
457     UInt_t pdgs[2];
458     
459     Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
460     pdgs[0]=211;
461     pdgs[1]=321;
462     Double_t minvD0 = d->InvMassD0();
463     pdgs[1]=211;
464     pdgs[0]=321;
465     Double_t minvD0bar = d->InvMassD0bar();
466     //cout<<"inside mass cut"<<endl;
467
468     Int_t pdgDgD0toKpi[2]={321,211};
469     Int_t lab=-9999;
470     if(fReadMC) lab=d->MatchToMC(421,mcArray,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
471     Double_t pt = d->Pt(); //mother pt
472  
473     if(lab>=0 && fReadMC){ //signal
474
475       //check pdg of the prongs
476       AliAODTrack *prong0=(AliAODTrack*)d->GetDaughter(0);
477       AliAODTrack *prong1=(AliAODTrack*)d->GetDaughter(1);
478       Int_t labprong[2];
479       labprong[0]=prong0->GetLabel();
480       labprong[1]=prong1->GetLabel();
481       AliAODMCParticle *mcprong=0;
482       Int_t PDGprong[2]={0,0};
483       for (Int_t iprong=0;iprong<2;iprong++){
484         if(labprong[iprong]>=0)  mcprong= (AliAODMCParticle*)mcArray->At(labprong[iprong]);
485         PDGprong[iprong]=mcprong->GetPdgCode();
486       }
487
488       //no mass cut ditributions: ptbis
489         
490       if(pt>0. && pt<=1.) {
491         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
492           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_1"))->Fill(d->PtProng(0));
493           ((TH1F*)fDistr->FindObject("hptKSnoMcut_1"))->Fill(d->PtProng(1));
494         }else {
495           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
496             ((TH1F*)fDistr->FindObject("hptKSnoMcut_1"))->Fill(d->PtProng(0));
497             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_1"))->Fill(d->PtProng(1));
498           }
499         }
500       }
501       if(pt>1. && pt<=2.) {
502         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
503           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_2"))->Fill(d->PtProng(0));
504           ((TH1F*)fDistr->FindObject("hptKSnoMcut_2"))->Fill(d->PtProng(1));
505         }else {
506           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
507             ((TH1F*)fDistr->FindObject("hptKSnoMcut_2"))->Fill(d->PtProng(0));
508             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_2"))->Fill(d->PtProng(1));
509           }
510         }
511       }
512       if(pt>2. && pt<=3.) {
513         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
514           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_3"))->Fill(d->PtProng(0));
515           ((TH1F*)fDistr->FindObject("hptKSnoMcut_3"))->Fill(d->PtProng(1));
516         }else {
517           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
518             ((TH1F*)fDistr->FindObject("hptKSnoMcut_3"))->Fill(d->PtProng(0));
519             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_3"))->Fill(d->PtProng(1));
520           }
521         }
522       }
523       if(pt>3. && pt<=5.) {
524
525         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
526           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_4"))->Fill(d->PtProng(0));
527           ((TH1F*)fDistr->FindObject("hptKSnoMcut_4"))->Fill(d->PtProng(1));
528         }else {
529           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
530             ((TH1F*)fDistr->FindObject("hptKSnoMcut_4"))->Fill(d->PtProng(0));
531             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_4"))->Fill(d->PtProng(1));
532           }
533         }
534       }
535       if(pt>5.)           {
536
537         if (TMath::Abs(PDGprong[0]) == 211 && TMath::Abs(PDGprong[1]) == 321){
538           ((TH1F*)fDistr->FindObject("hptpiSnoMcut_5"))->Fill(d->PtProng(0));
539           ((TH1F*)fDistr->FindObject("hptKSnoMcut_5"))->Fill(d->PtProng(1));
540         }else {
541           if (TMath::Abs(PDGprong[0]) == 321 && TMath::Abs(PDGprong[1]) == 211){
542             ((TH1F*)fDistr->FindObject("hptKSnoMcut_5"))->Fill(d->PtProng(0));
543             ((TH1F*)fDistr->FindObject("hptpiSnoMcut_5"))->Fill(d->PtProng(1));
544           }
545         }
546       }
547
548       if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421){//D0
549                          
550         //no mass cut ditributions: mass, costhetastar
551         
552         if(pt>0. && pt<=1.) {
553           ((TH1F*)fDistr->FindObject("hMassS_1"))->Fill(minvD0);
554         }
555         if(pt>1. && pt<=2.) {
556           ((TH1F*)fDistr->FindObject("hMassS_2"))->Fill(minvD0);
557
558         }
559         if(pt>2. && pt<=3.) {
560           ((TH1F*)fDistr->FindObject("hMassS_3"))->Fill(minvD0);
561
562         }
563         if(pt>3. && pt<=5.) {
564           ((TH1F*)fDistr->FindObject("hMassS_4"))->Fill(minvD0);
565
566         }
567         if(pt>5.)           {
568           ((TH1F*)fDistr->FindObject("hMassS_5"))->Fill(minvD0);
569
570         }
571
572       }
573       else { //D0bar
574
575         //no mass cut ditributions: mass
576
577         if(pt>0. && pt<=1.) {
578           ((TH1F*)fDistr->FindObject("hMassS_1"))->Fill(minvD0bar);
579
580         }
581         if(pt>1. && pt<=2.) {
582           ((TH1F*)fDistr->FindObject("hMassS_2"))->Fill(minvD0bar);
583
584         }
585         if(pt>2. && pt<=3.) {
586           ((TH1F*)fDistr->FindObject("hMassS_3"))->Fill(minvD0bar);
587
588         }
589         if(pt>3. && pt<=5.) {
590           ((TH1F*)fDistr->FindObject("hMassS_4"))->Fill(minvD0bar);
591
592         }
593         if(pt>5.)           {
594           ((TH1F*)fDistr->FindObject("hMassS_5"))->Fill(minvD0bar);
595         }
596
597       }
598
599       //apply cut on invariant mass on the pair
600       if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
601
602         if(fArray==1) cout<<"LS signal: ERROR"<<endl;
603         for (Int_t iprong=0; iprong<2; iprong++){
604           AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(iprong);
605           Int_t labprong=prong->GetLabel();
606
607           //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
608           AliAODMCParticle *mcprong=0;
609           if(labprong>=0)  mcprong= (AliAODMCParticle*)mcArray->At(labprong);
610           Int_t pdgprong=mcprong->GetPdgCode();
611           if(TMath::Abs(pdgprong)==211) {
612             //cout<<"pi"<<endl;
613             if(pt>0. && pt<=1.){
614               ((TH1F*)fDistr->FindObject("hptpiS_1"))->Fill(d->PtProng(iprong));
615               ((TH1F*)fDistr->FindObject("hd0piS_1"))->Fill(d->Getd0Prong(iprong));
616             }
617             if(pt>1. && pt<=2.){
618               ((TH1F*)fDistr->FindObject("hptpiS_2"))->Fill(d->PtProng(iprong));
619               ((TH1F*)fDistr->FindObject("hd0piS_2"))->Fill(d->Getd0Prong(iprong));
620             }
621
622             if(pt>2. && pt<=3.){
623               ((TH1F*)fDistr->FindObject("hptpiS_3"))->Fill(d->PtProng(iprong));
624               ((TH1F*)fDistr->FindObject("hd0piS_3"))->Fill(d->Getd0Prong(iprong));
625
626             }
627             if(pt>3. && pt<=5.){
628               ((TH1F*)fDistr->FindObject("hptpiS_4"))->Fill(d->PtProng(iprong));
629               ((TH1F*)fDistr->FindObject("hd0piS_4"))->Fill(d->Getd0Prong(iprong));
630
631             }
632             if(pt>5.)          {
633               ((TH1F*)fDistr->FindObject("hptpiS_5"))->Fill(d->PtProng(iprong));
634               ((TH1F*)fDistr->FindObject("hd0piS_5"))->Fill(d->Getd0Prong(iprong));
635
636             }
637
638
639         }
640
641         if(TMath::Abs(pdgprong)==321) {
642             //cout<<"kappa"<<endl;
643             if(pt>0. && pt<=1.){
644               ((TH1F*)fDistr->FindObject("hptKS_1"))->Fill(d->PtProng(iprong));
645               ((TH1F*)fDistr->FindObject("hd0KS_1"))->Fill(d->Getd0Prong(iprong));
646             }
647             if(pt>1. && pt<=2.){
648               ((TH1F*)fDistr->FindObject("hptKS_2"))->Fill(d->PtProng(iprong));
649               ((TH1F*)fDistr->FindObject("hd0KS_2"))->Fill(d->Getd0Prong(iprong));
650             }
651             if(pt>2. && pt<=3.){
652               ((TH1F*)fDistr->FindObject("hptKS_3"))->Fill(d->PtProng(iprong));
653               ((TH1F*)fDistr->FindObject("hd0KS_3"))->Fill(d->Getd0Prong(iprong));
654             }
655             if(pt>3. && pt<=5.){
656               ((TH1F*)fDistr->FindObject("hptKS_4"))->Fill(d->PtProng(iprong));
657               ((TH1F*)fDistr->FindObject("hd0KS_4"))->Fill(d->Getd0Prong(iprong));
658
659             }
660             if(pt>5.)          {
661               ((TH1F*)fDistr->FindObject("hptKS_5"))->Fill(d->PtProng(iprong));
662               ((TH1F*)fDistr->FindObject("hd0KS_5"))->Fill(d->Getd0Prong(iprong));
663
664             }
665
666           }
667
668           if(pt>0. && pt<=1.){
669             ((TH1F*)fDistr->FindObject("hdcaS_1"))->Fill(d->GetDCA());
670
671             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
672               ((TH1F*)fDistr->FindObject("hcosthetastarS_1"))->Fill(d->CosThetaStarD0());
673             else ((TH1F*)fDistr->FindObject("hcosthetastarS_1"))->Fill(d->CosThetaStarD0bar());
674
675           }
676           if(pt>1. && pt<=2.){
677             ((TH1F*)fDistr->FindObject("hdcaS_2"))->Fill(d->GetDCA());
678
679             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
680               ((TH1F*)fDistr->FindObject("hcosthetastarS_2"))->Fill(d->CosThetaStarD0());
681             else ((TH1F*)fDistr->FindObject("hcosthetastarS_2"))->Fill(d->CosThetaStarD0bar());
682
683           }
684           if(pt>2. && pt<=3.){
685             ((TH1F*)fDistr->FindObject("hdcaS_3"))->Fill(d->GetDCA());
686
687             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
688               ((TH1F*)fDistr->FindObject("hcosthetastarS_3"))->Fill(d->CosThetaStarD0());
689             else ((TH1F*)fDistr->FindObject("hcosthetastarS_3"))->Fill(d->CosThetaStarD0bar());
690
691           }
692           if(pt>3. && pt<=5.){
693             ((TH1F*)fDistr->FindObject("hdcaS_4"))->Fill(d->GetDCA());
694
695             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
696               ((TH1F*)fDistr->FindObject("hcosthetastarS_4"))->Fill(d->CosThetaStarD0());
697             else ((TH1F*)fDistr->FindObject("hcosthetastarS_4"))->Fill(d->CosThetaStarD0bar());
698
699           }
700           if(pt>5.)          {
701             ((TH1F*)fDistr->FindObject("hdcaS_5"))->Fill(d->GetDCA());
702
703             if (((AliAODMCParticle*)mcArray->At(lab))->GetPdgCode() == 421)
704               ((TH1F*)fDistr->FindObject("hcosthetastarS_5"))->Fill(d->CosThetaStarD0());
705             else ((TH1F*)fDistr->FindObject("hcosthetastarS_5"))->Fill(d->CosThetaStarD0bar());
706
707           }
708           
709         }
710
711         if(pt>0. && pt<=1.){
712           ((TH1F*)fDistr->FindObject("hd0d0S_1"))->Fill(d->Prodd0d0());
713           ((TH1F*)fDistr->FindObject("hcosthetapointS_1"))->Fill(d->CosPointingAngle());
714           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_1"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
715         }
716         if(pt>1. && pt<=2.){
717           ((TH1F*)fDistr->FindObject("hd0d0S_2"))->Fill(d->Prodd0d0());
718           ((TH1F*)fDistr->FindObject("hcosthetapointS_2"))->Fill(d->CosPointingAngle());
719           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_2"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
720         }
721         if(pt>2. && pt<=3.){
722           ((TH1F*)fDistr->FindObject("hd0d0S_3"))->Fill(d->Prodd0d0());
723           ((TH1F*)fDistr->FindObject("hcosthetapointS_3"))->Fill(d->CosPointingAngle());
724           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_3"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
725         }
726         if(pt>3. && pt<=5.){
727           ((TH1F*)fDistr->FindObject("hd0d0S_4"))->Fill(d->Prodd0d0());
728           ((TH1F*)fDistr->FindObject("hcosthetapointS_4"))->Fill(d->CosPointingAngle());
729           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_4"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
730         }
731         if(pt>5.)          {
732           ((TH1F*)fDistr->FindObject("hd0d0S_5"))->Fill(d->Prodd0d0());
733           ((TH1F*)fDistr->FindObject("hcosthetapointS_5"))->Fill(d->CosPointingAngle());
734           ((TH1F*)fDistr->FindObject("hcosthpointd0d0S_5"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
735         }
736
737       } //invmass cut
738
739     } else{ //Background or LS
740       //cout<<"is background"<<endl;
741       Double_t pt = d->Pt();
742
743       //no mass cut distributions: mass, ptbis
744       if(pt>0. && pt<=1.) {
745         ((TH1F*)fDistr->FindObject("hMassB_1"))->Fill(minvD0);
746         ((TH1F*)fDistr->FindObject("hMassB_1"))->Fill(minvD0bar);
747
748         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_1"))->Fill(d->PtProng(0));
749         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_1"))->Fill(d->PtProng(0));
750         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_1"))->Fill(d->PtProng(1));
751
752       }
753       if(pt>1. && pt<=2.) {
754         ((TH1F*)fDistr->FindObject("hMassB_2"))->Fill(minvD0);
755         ((TH1F*)fDistr->FindObject("hMassB_2"))->Fill(minvD0bar);
756
757         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_2"))->Fill(d->PtProng(0));
758         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_2"))->Fill(d->PtProng(0));
759         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_2"))->Fill(d->PtProng(1));
760
761       }
762       if(pt>2. && pt<=3.) {
763         ((TH1F*)fDistr->FindObject("hMassB_3"))->Fill(minvD0);
764         ((TH1F*)fDistr->FindObject("hMassB_3"))->Fill(minvD0bar);
765
766         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_3"))->Fill(d->PtProng(0));
767         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_3"))->Fill(d->PtProng(0));
768         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_3"))->Fill(d->PtProng(1));
769
770       }
771       if(pt>3. && pt<=5.) {
772         ((TH1F*)fDistr->FindObject("hMassB_4"))->Fill(minvD0);
773         ((TH1F*)fDistr->FindObject("hMassB_4"))->Fill(minvD0bar);
774
775         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_4"))->Fill(d->PtProng(0));
776         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_4"))->Fill(d->PtProng(0));
777         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_4"))->Fill(d->PtProng(1));
778  
779       }
780       if(pt>5.)           { 
781         ((TH1F*)fDistr->FindObject("hMassB_5"))->Fill(minvD0);
782         ((TH1F*)fDistr->FindObject("hMassB_5"))->Fill(minvD0bar);
783
784         ((TH1F*)fDistr->FindObject("hptB1prongnoMcut_5"))->Fill(d->PtProng(0));
785         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_5"))->Fill(d->PtProng(0));
786         ((TH1F*)fDistr->FindObject("hptB2prongsnoMcut_5"))->Fill(d->PtProng(1));
787  
788       }
789
790       //apply cut on invariant mass on the pair
791       if(TMath::Abs(minvD0-mPDG)<0.03 || TMath::Abs(minvD0bar-mPDG)<0.03){
792
793
794         AliAODTrack *prong=(AliAODTrack*)d->GetDaughter(0);
795         if(!prong) cout<<"No daughter found";
796         else{
797           if(prong->Charge()==1) {fTotPosPairs[4]++;} else {fTotNegPairs[4]++;}
798         }
799         if(pt>0. && pt<=1.){
800           //normalise pt distr to half afterwards
801           ((TH1F*)fDistr->FindObject("hptB_1"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_1"))->Fill(d->PtProng(1));
802           ((TH1F*)fDistr->FindObject("hd0B_1"))->Fill(d->Getd0Prong(0));
803           ((TH1F*)fDistr->FindObject("hdcaB_1"))->Fill(d->GetDCA());
804           ((TH1F*)fDistr->FindObject("hcosthetastarB_1"))->Fill(d->CosThetaStarD0());
805           ((TH1F*)fDistr->FindObject("hcosthetastarB_1"))->Fill(d->CosThetaStarD0bar());        
806           ((TH1F*)fDistr->FindObject("hd0d0B_1"))->Fill(d->Prodd0d0());
807           ((TH1F*)fDistr->FindObject("hcosthetapointB_1"))->Fill(d->CosPointingAngle());
808           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_1"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
809         }
810         if(pt>1. && pt<=2.){
811           ((TH1F*)fDistr->FindObject("hptB_2"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_2"))->Fill(d->PtProng(1));
812           ((TH1F*)fDistr->FindObject("hd0B_2"))->Fill(d->Getd0Prong(0));
813           ((TH1F*)fDistr->FindObject("hdcaB_2"))->Fill(d->GetDCA());
814           ((TH1F*)fDistr->FindObject("hcosthetastarB_2"))->Fill(d->CosThetaStarD0());
815           ((TH1F*)fDistr->FindObject("hcosthetastarB_2"))->Fill(d->CosThetaStarD0bar());        
816           ((TH1F*)fDistr->FindObject("hd0d0B_2"))->Fill(d->Prodd0d0());
817           ((TH1F*)fDistr->FindObject("hcosthetapointB_2"))->Fill(d->CosPointingAngle());
818           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_2"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
819         }
820         if(pt>2. && pt<=3.){
821           ((TH1F*)fDistr->FindObject("hptB_3"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_3"))->Fill(d->PtProng(1));
822           ((TH1F*)fDistr->FindObject("hd0B_3"))->Fill(d->Getd0Prong(0));
823           ((TH1F*)fDistr->FindObject("hdcaB_3"))->Fill(d->GetDCA());
824           ((TH1F*)fDistr->FindObject("hcosthetastarB_3"))->Fill(d->CosThetaStarD0());
825           ((TH1F*)fDistr->FindObject("hcosthetastarB_3"))->Fill(d->CosThetaStarD0bar());        
826           ((TH1F*)fDistr->FindObject("hd0d0B_3"))->Fill(d->Prodd0d0());
827           ((TH1F*)fDistr->FindObject("hcosthetapointB_3"))->Fill(d->CosPointingAngle());
828           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_3"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
829         }
830         if(pt>3. && pt<=5.){
831           ((TH1F*)fDistr->FindObject("hptB_4"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_4"))->Fill(d->PtProng(1));
832           ((TH1F*)fDistr->FindObject("hd0B_4"))->Fill(d->Getd0Prong(0));
833           ((TH1F*)fDistr->FindObject("hdcaB_4"))->Fill(d->GetDCA());
834           ((TH1F*)fDistr->FindObject("hcosthetastarB_4"))->Fill(d->CosThetaStarD0());
835           ((TH1F*)fDistr->FindObject("hcosthetastarB_4"))->Fill(d->CosThetaStarD0bar());        
836           ((TH1F*)fDistr->FindObject("hd0d0B_4"))->Fill(d->Prodd0d0());
837           ((TH1F*)fDistr->FindObject("hcosthetapointB_4"))->Fill(d->CosPointingAngle());
838           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_4"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
839         }
840         if(pt>5.)         {
841           ((TH1F*)fDistr->FindObject("hptB_5"))->Fill(d->PtProng(0));((TH1F*)fDistr->FindObject("hptB_5"))->Fill(d->PtProng(1));
842           ((TH1F*)fDistr->FindObject("hd0B_5"))->Fill(d->Getd0Prong(0));
843           ((TH1F*)fDistr->FindObject("hdcaB_5"))->Fill(d->GetDCA());
844           ((TH1F*)fDistr->FindObject("hcosthetastarB_5"))->Fill(d->CosThetaStarD0());
845           ((TH1F*)fDistr->FindObject("hcosthetastarB_5"))->Fill(d->CosThetaStarD0bar());        
846           ((TH1F*)fDistr->FindObject("hd0d0B_5"))->Fill(d->Prodd0d0());
847           ((TH1F*)fDistr->FindObject("hcosthetapointB_5"))->Fill(d->CosPointingAngle());
848           ((TH1F*)fDistr->FindObject("hcosthpointd0d0B_5"))->Fill(d->CosPointingAngle(),d->Prodd0d0());
849         }       
850
851
852       }// end if inv mass cut
853     }//end if background
854
855   
856   
857      //cuts order
858 //       printf("    |M-MD0| [GeV]    < %f\n",fD0toKpiCuts[0]);
859 //     printf("    dca    [cm]  < %f\n",fD0toKpiCuts[1]);
860 //     printf("    cosThetaStar     < %f\n",fD0toKpiCuts[2]);
861 //     printf("    pTK     [GeV/c]    > %f\n",fD0toKpiCuts[3]);
862 //     printf("    pTpi    [GeV/c]    > %f\n",fD0toKpiCuts[4]);
863 //     printf("    |d0K|  [cm]  < %f\n",fD0toKpiCuts[5]);
864 //     printf("    |d0pi| [cm]  < %f\n",fD0toKpiCuts[6]);
865 //     printf("    d0d0  [cm^2] < %f\n",fD0toKpiCuts[7]);
866 //     printf("    cosThetaPoint    > %f\n",fD0toKpiCuts[8]);
867   
868     Int_t ptbin=0;
869     Bool_t isInRange=kFALSE;
870     
871     //cout<<"P_t = "<<pt<<endl;
872     if (pt>0. && pt<=1.) {
873       ptbin=1;
874       isInRange=kTRUE;
875       /*
876       //test d0 cut
877       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.1,-0.0002,0.5);
878       fVHFmycuts->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.1,-0.00025,0.7);
879       */
880       //original
881       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.5);
882       fVHFmycuts->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.00025,0.7);
883       
884 //       fVHFPPR->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,0.05,0.05,-0.0002,0.7);
885 //       fVHFmycuts->SetD0toKpiCuts(0.7,0.04,0.8,0.5,0.5,1,1,-0.00015,0.5);
886       //printf("I'm in the bin %d\n",ptbin);
887       
888     }
889     
890     if(pt>1. && pt<=3.) {
891       if(pt>1. && pt<=2.) ptbin=2;  
892       if(pt>2. && pt<=3.) ptbin=3;  
893       isInRange=kTRUE;
894       /*
895       //test d0 cut
896       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.0002,0.6);
897       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,0.1,-0.00025,0.8);
898       */
899       //original
900       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0002,0.6);
901       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,1,1,-0.00025,0.8);
902       
903       //printf("I'm in the bin %d\n",ptbin);
904     }
905  
906     if(pt>3. && pt<=5.){
907       ptbin=4;  
908       isInRange=kTRUE;
909       /*
910       //test d0 cut
911       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.0001,0.8);
912       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.00015,0.8);
913       */
914       //original
915       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.0001,0.8);
916       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.8);
917       
918       //printf("I'm in the bin %d\n",ptbin);
919     }
920     if(pt>5.&& pt<=10.){ //additional upper limit to compare with Correction Framework
921       ptbin=5;
922       isInRange=kTRUE;
923       /*
924       //test d0 cut
925       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.00005,0.8);
926       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.1,-0.00015,0.9);
927       */
928       //original
929       fVHFPPR->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00005,0.8);
930       fVHFmycuts->SetD0toKpiCuts(0.7,0.02,0.8,0.7,0.7,0.05,0.05,-0.00015,0.9);
931       
932     }//if(pt>5)    if (pt>0. && pt<=1.) {
933     //printf("I'm in the bin %d\n",ptbin);
934     //old
935     //fVHF->SetD0toKpiCuts(0.7,0.03,0.8,0.06,0.06,0.05,0.05,-0.0002,0.6); //2.p-p vertex reconstructed    
936     if(prongsinacc && isInRange){
937       FillHists(ptbin,d,mcArray,fVHFPPR,fOutputPPR);
938       FillHists(ptbin,d,mcArray,fVHFmycuts,fOutputmycuts);
939     }
940     if(unsetvtx) d->UnsetOwnPrimaryVtx();
941   } //end for prongs
942   
943      
944   
945    
946   // Post the data
947   PostData(1,fOutputPPR);
948   PostData(2,fOutputmycuts);
949   PostData(4,fDistr);
950   return;
951 }
952 //____________________________________________________________________________*
953 void AliAnalysisTaskSED0Mass::FillHists(Int_t ptbin, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAnalysisVertexingHF *vhf, TList *listout){
954   //
955   // function used in UserExec:
956   //
957
958
959   Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
960
961   Int_t okD0=0,okD0bar=0;
962   //cout<<"inside FillHist"<<endl;
963   if(part->SelectD0(vhf->GetD0toKpiCuts(),okD0,okD0bar)) {//selected
964     Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
965     //printf("SELECTED\n");
966
967
968     AliAODTrack *prong=(AliAODTrack*)part->GetDaughter(0);
969     if(!prong) cout<<"No daughter found";
970     else{
971       if(prong->Charge()==1) {fTotPosPairs[ptbin-1]++;} else {fTotNegPairs[ptbin-1]++;}
972     }
973
974     TString fillthis="";
975     Int_t pdgDgD0toKpi[2]={321,211};
976     Int_t labD0=-1;
977     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)
978
979     //count candidates selected by cuts
980     fNentries->Fill(2);
981     //count true D0 selected by cuts
982     if (fReadMC && labD0>=0) fNentries->Fill(3);
983     PostData(3,fNentries);
984
985     if (okD0==1) {
986       fillthis="histMass_";
987       fillthis+=ptbin;
988       //cout<<"Filling "<<fillthis<<endl;
989
990       //printf("Fill mass with D0");
991       ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
992  
993       if(labD0>=0) {
994         if(fArray==1) cout<<"LS signal ERROR"<<endl;
995
996         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
997         Int_t pdgD0 = partD0->GetPdgCode();
998         //cout<<"pdg = "<<pdgD0<<endl;
999         if (pdgD0==421){ //D0
1000           //cout<<"Fill S with D0"<<endl;
1001           fillthis="histSgn_";
1002           fillthis+=ptbin;
1003           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1004           if(TMath::Abs(invmassD0 - mPDG) < 0.027){
1005             fillthis="histSgn27_";
1006             fillthis+=ptbin;
1007             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1008           }
1009         } else{ //it was a D0bar
1010           fillthis="histRfl_";
1011           fillthis+=ptbin;
1012           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1013         }
1014       } else {//background
1015         fillthis="histBkg_";
1016         fillthis+=ptbin;
1017         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0);
1018       }
1019       
1020     }
1021     if (okD0bar==1) {
1022       fillthis="histMass_";
1023       fillthis+=ptbin;
1024       //printf("Fill mass with D0bar");
1025       ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar);
1026  
1027       if(labD0>=0) {
1028         if(fArray==1) cout<<"LS signal ERROR"<<endl;
1029         AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1030         Int_t pdgD0 = partD0->GetPdgCode();
1031         //cout<<" pdg = "<<pdgD0<<endl;
1032         if (pdgD0==-421){ //D0bar
1033           fillthis="histSgn_";
1034           fillthis+=ptbin;
1035           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1036           if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
1037             fillthis="histSgn27_";
1038             fillthis+=ptbin;
1039             ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1040           }
1041
1042           
1043         } else{
1044           fillthis="histRfl_";
1045           fillthis+=ptbin;
1046           ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1047         }
1048       } else {//background or LS
1049         fillthis="histBkg_";
1050         fillthis+=ptbin;
1051         ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
1052       }
1053     }
1054
1055
1056   } //else cout<<"NOT SELECTED"<<endl;
1057
1058 }
1059 //________________________________________________________________________
1060 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
1061 {
1062   // Terminate analysis
1063   //
1064   if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
1065
1066
1067   fOutputPPR = dynamic_cast<TList*> (GetOutputData(1));
1068   if (!fOutputPPR) {     
1069     printf("ERROR: fOutputthight not available\n");
1070     return;
1071   }
1072   fOutputmycuts = dynamic_cast<TList*> (GetOutputData(2));
1073   if (!fOutputmycuts) {     
1074     printf("ERROR: fOutputmycuts not available\n");
1075     return;
1076   }
1077   fDistr = dynamic_cast<TList*> (GetOutputData(4));
1078   if (!fDistr) {
1079     printf("ERROR: fDistr not available\n");
1080     return;
1081   }
1082
1083   if(fArray==1){
1084     for(Int_t ipt=0;ipt<4;ipt++){
1085       fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[ipt]*fTotNegPairs[ipt]); 
1086
1087
1088       if(fLsNormalization>0) {
1089         
1090         TString massName="histMass_";
1091         massName+=ipt+1;
1092         ((TH1F*)fOutputPPR->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputPPR->FindObject(massName))->GetEntries());
1093         ((TH1F*)fOutputmycuts->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputmycuts->FindObject(massName))->GetEntries());
1094       }
1095     }
1096
1097     fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
1098
1099     if(fLsNormalization>0) {
1100
1101       TString nameDistr="hptB";
1102       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1103       nameDistr="hdcaB";
1104       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1105       nameDistr="hcosthetastarB";
1106       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1107       nameDistr="hd0B";
1108       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1109       nameDistr="hd0d0B";
1110       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1111       nameDistr="hcosthetapointB";
1112       ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
1113       nameDistr="hcosthpointd0d0B";
1114       ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
1115
1116     }
1117   }
1118
1119   TString cvname;
1120
1121   if (fArray==0){
1122     cvname="D0invmass";
1123   } else cvname="LSinvmass";
1124
1125   TCanvas *c1=new TCanvas(cvname,cvname);
1126   c1->cd();
1127   ((TH1F*)fOutputPPR->FindObject("histMass_4"))->Draw();
1128
1129
1130   return;
1131 }
1132