Changes to compile with Root6 on macosx64
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEDplus.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //*************************************************************************
19 // Class AliAnalysisTaskSEDplus
20 // AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and 
21 //comparison of heavy-flavour decay candidates
22 // to MC truth (kinematics stored in the AOD)
23 // Authors: Renu Bala, bala@to.infn.it
24 // F. Prino, prino@to.infn.it
25 // G. Ortona, ortona@to.infn.it
26 /////////////////////////////////////////////////////////////
27
28 #include <TClonesArray.h>
29 #include <TNtuple.h>
30 #include <TCanvas.h>
31 #include <TList.h>
32 #include <TString.h>
33 #include <TDatabasePDG.h>
34
35 #include "AliAnalysisManager.h"
36 #include "AliRDHFCutsDplustoKpipi.h"
37 #include "AliAODHandler.h"
38 #include "AliAODEvent.h"
39 #include "AliAODVertex.h"
40 #include "AliAODTrack.h"
41 #include "AliAODRecoDecayHF3Prong.h"
42 #include "AliAnalysisVertexingHF.h"
43 #include "AliAnalysisTaskSE.h"
44 #include "AliAnalysisTaskSEDplus.h"
45 #include "AliNormalizationCounter.h"
46 #include "AliVertexingHFUtils.h"
47 ClassImp(AliAnalysisTaskSEDplus)
48
49 //________________________________________________________________________
50 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
51 AliAnalysisTaskSE(),
52   fOutput(0), 
53   fHistNEvents(0),
54   fPtVsMass(0),
55   fPtVsMassTC(0),
56   fYVsPt(0),
57   fYVsPtTC(0),
58   fYVsPtSig(0),
59   fYVsPtSigTC(0),
60   fPhiEtaCand(0),
61   fPhiEtaCandSigReg(0),
62   fSPDMult(0),
63   fNtupleDplus(0),
64   fUpmasslimit(1.965),
65   fLowmasslimit(1.765),
66   fNPtBins(0),
67   fBinWidth(0.002),
68   fListCuts(0),
69   fRDCutsAnalysis(0),
70   fCounter(0),
71   fFillNtuple(kFALSE),
72   fReadMC(kFALSE),
73   fUseStrangeness(kFALSE),
74   fUseBit(kTRUE),
75   fCutsDistr(kFALSE),
76   fDoImpPar(kFALSE),
77   fNImpParBins(400),
78   fLowerImpPar(-1000.),
79   fHigherImpPar(1000.),
80   fDoLS(0),
81   fEtaSelection(0),
82   fSystem(0)
83 {
84   // Default constructor
85
86   for(Int_t i=0;i<3;i++){
87     if(fHistCentrality[i])fHistCentrality[i]=0;
88     if(fCorreld0Kd0pi[i])fCorreld0Kd0pi[i]=0;
89   }
90   
91   for(Int_t i=0; i<5; i++)fHistMassPtImpParTC[i]=0;
92     
93   
94   for(Int_t i=0;i<3*kMaxPtBins;i++){
95     if(fMassHist[i])fMassHist[i]=0;
96     if(fCosPHist[i])fCosPHist[i]=0;
97     if(fDLenHist[i])fDLenHist[i]=0;
98     if(fSumd02Hist[i])fSumd02Hist[i]=0;
99     if(fSigVertHist[i])fSigVertHist[i]=0;
100     if(fPtMaxHist[i])fPtMaxHist[i]=0;
101     if(fPtKHist[i])fPtKHist[i]=0;
102     if(fPtpi1Hist[i])fPtpi1Hist[i]=0;
103     if(fPtpi2Hist[i])fPtpi2Hist[i]=0;
104     if(fDCAHist[i])fDCAHist[i]=0;
105     if(fMassHistTC[i])fMassHistTC[i]=0;
106     if(fMassHistTCPlus[i])fMassHistTCPlus[i]=0;
107     if(fMassHistTCMinus[i])fMassHistTCMinus[i]=0;
108     
109     if(fDLxy[i])fDLxy[i]=0;
110     if(fDLxyTC[i])fDLxyTC[i]=0;
111     if(fCosxy[i])fCosxy[i]=0;
112     if(fCosxyTC[i])fCosxyTC[i]=0;
113     if(fCosPHistLS[i])fCosPHistLS[i]=0;
114     if(fDLenHistLS[i])fDLenHistLS[i]=0;
115     if(fSumd02HistLS[i])fSumd02HistLS[i]=0;
116     if(fSigVertHistLS[i])fSigVertHistLS[i]=0;
117     if(fPtMaxHistLS[i])fPtMaxHistLS[i]=0;
118     if(fDCAHistLS[i])fDCAHistLS[i]=0;
119   }
120   for(Int_t i=0;i<5*kMaxPtBins;i++){
121     if(fMassHistLS[i])fMassHistLS[i]=0;
122     if(fMassHistLSTC[i])fMassHistLSTC[i]=0;
123   }
124   for(Int_t i=0;i<kMaxPtBins+1;i++){
125     fArrayBinLimits[i]=0;
126   }
127
128 }
129
130 //________________________________________________________________________
131 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,Bool_t fillNtuple):
132   AliAnalysisTaskSE(name),
133   fOutput(0),
134   fHistNEvents(0),
135   fPtVsMass(0),
136   fPtVsMassTC(0),
137   fYVsPt(0),
138   fYVsPtTC(0),
139   fYVsPtSig(0),
140   fYVsPtSigTC(0),
141   fPhiEtaCand(0),
142   fPhiEtaCandSigReg(0),
143   fSPDMult(0),
144   fNtupleDplus(0),
145   fUpmasslimit(1.965),
146   fLowmasslimit(1.765),
147   fNPtBins(0),
148   fBinWidth(0.002),
149   fListCuts(0),
150   fRDCutsAnalysis(dpluscutsana),
151   fCounter(0),
152   fFillNtuple(fillNtuple),
153   fReadMC(kFALSE),
154   fUseStrangeness(kFALSE),
155   fUseBit(kTRUE),
156   fCutsDistr(kFALSE),
157   fDoImpPar(kFALSE),
158   fNImpParBins(400),
159   fLowerImpPar(-1000.),
160   fHigherImpPar(1000.),
161   fDoLS(0),
162   fEtaSelection(0),
163   fSystem(0)
164 {
165   // 
166   // Standrd constructor
167   //
168   fNPtBins=fRDCutsAnalysis->GetNPtBins();
169   
170   for(Int_t i=0;i<3;i++){
171     if(fHistCentrality[i])fHistCentrality[i]=0;
172     if(fCorreld0Kd0pi[i])fCorreld0Kd0pi[i]=0;
173   }
174   
175   for(Int_t i=0; i<5; i++)fHistMassPtImpParTC[i]=0;
176     
177   for(Int_t i=0;i<3*kMaxPtBins;i++){
178     if(fMassHist[i])fMassHist[i]=0;
179     if(fCosPHist[i])fCosPHist[i]=0;
180     if(fDLenHist[i])fDLenHist[i]=0;
181     if(fSumd02Hist[i])fSumd02Hist[i]=0;
182     if(fSigVertHist[i])fSigVertHist[i]=0;
183     if(fPtMaxHist[i])fPtMaxHist[i]=0;
184     if(fPtKHist[i])fPtKHist[i]=0;
185     if(fPtpi1Hist[i])fPtpi1Hist[i]=0;
186     if(fPtpi2Hist[i])fPtpi2Hist[i]=0;
187     if(fDCAHist[i])fDCAHist[i]=0;
188     if(fMassHistTC[i])fMassHistTC[i]=0;
189     if(fMassHistTCPlus[i])fMassHistTCPlus[i]=0;
190     if(fMassHistTCMinus[i])fMassHistTCMinus[i]=0;
191     
192     if(fDLxy[i])fDLxy[i]=0;
193     if(fDLxyTC[i])fDLxyTC[i]=0;
194     if(fCosxy[i])fCosxy[i]=0;
195     if(fCosxyTC[i])fCosxyTC[i]=0;
196     if(fCosPHistLS[i])fCosPHistLS[i]=0;
197     if(fDLenHistLS[i])fDLenHistLS[i]=0;
198     if(fSumd02HistLS[i])fSumd02HistLS[i]=0;
199     if(fSigVertHistLS[i])fSigVertHistLS[i]=0;
200     if(fPtMaxHistLS[i])fPtMaxHistLS[i]=0;
201     if(fDCAHistLS[i])fDCAHistLS[i]=0;
202   }
203   for(Int_t i=0;i<5*kMaxPtBins;i++){
204     if(fMassHistLS[i])fMassHistLS[i]=0;
205     if(fMassHistLSTC[i])fMassHistLSTC[i]=0;
206   }
207   for(Int_t i=0;i<kMaxPtBins+1;i++){
208     fArrayBinLimits[i]=0;
209   }
210   
211   
212   // Default constructor
213   // Output slot #1 writes into a TList container
214   DefineOutput(1,TList::Class());  //My private output
215   // Output slot #2 writes cut to private output
216   //  DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
217   DefineOutput(2,TList::Class());
218   // Output slot #3 writes cut to private output
219   DefineOutput(3,AliNormalizationCounter::Class());
220   
221   if(fFillNtuple){
222     // Output slot #4 writes into a TNtuple container
223     DefineOutput(4,TNtuple::Class());  //My private output
224   }
225 }
226
227 //________________________________________________________________________
228 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
229 {
230   //
231   // Destructor
232   //
233   if (fOutput) {
234     delete fOutput;
235     fOutput = 0;
236   }
237
238   
239   if(fNtupleDplus){
240     delete fNtupleDplus;
241     fNtupleDplus=0;
242   }
243   if (fListCuts) {
244     delete fListCuts;
245     fListCuts = 0;
246   }
247   if(fRDCutsAnalysis){
248     delete fRDCutsAnalysis;
249     fRDCutsAnalysis = 0;
250   }
251   if(fCounter){
252     delete fCounter;
253     fCounter = 0;
254   }
255
256
257 }  
258 //_________________________________________________________________
259 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
260   // set invariant mass limits
261   Float_t bw=GetBinWidth();
262   fUpmasslimit = 1.865+range;
263   fLowmasslimit = 1.865-range;
264   SetBinWidth(bw);
265 }
266 //_________________________________________________________________
267 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
268   // set invariant mass limits
269   if(uplimit>lowlimit)
270     {
271       Float_t bw=GetBinWidth();
272       fUpmasslimit = lowlimit;
273       fLowmasslimit = uplimit;
274       SetBinWidth(bw);
275     }
276 }
277 //________________________________________________________________
278 void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
279   Float_t width=w;
280   Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
281   Int_t missingbins=4-nbins%4;
282   nbins=nbins+missingbins;
283   width=(fUpmasslimit-fLowmasslimit)/nbins;
284   if(missingbins!=0){
285     printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
286   }
287   else{
288     if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
289   }
290   fBinWidth=width;
291 }
292 //_________________________________________________________________
293 Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
294   return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
295 }
296 //_________________________________________________________________
297 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
298   //
299   //
300   // Fill the Like Sign histograms
301   //
302   if(fDebug>1)printf("started LS\n");
303   
304   //histograms for like sign
305   Int_t nbins=GetNBinsHistos();;
306   TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
307   TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
308   
309   Int_t nPosTrks=0,nNegTrks=0;
310   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
311   Int_t nDplusLS=0;
312   Int_t nDminusLS=0;
313   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
314   Int_t index=0; 
315   
316   // loop over like sign candidates
317   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
318     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
319     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
320     Bool_t unsetvtx=kFALSE;
321     if(!d->GetOwnPrimaryVtx()) {
322       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
323       unsetvtx=kTRUE;
324     }
325     Double_t ptCand = d->Pt();
326     Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
327     if(iPtBin<0)continue;
328     Int_t sign= d->GetCharge();
329     if(sign>0){
330       nPosTrks++;
331     }else{
332       nNegTrks++;
333     }
334     fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
335     Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
336     if(passTightCuts>0){
337       Float_t invMass = d->InvMassDplus();
338       index=GetLSHistoIndex(iPtBin);
339       fMassHistLS[index+1]->Fill(invMass);
340       if(sign>0){
341         histLSPlus->Fill(invMass);
342         nDplusLS++;
343       }else{
344         histLSMinus->Fill(invMass);
345         nDminusLS++;
346       }
347       if(fCutsDistr){
348         Double_t dlen=d->DecayLength();
349         Double_t cosp=d->CosPointingAngle();
350         Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
351         Double_t dca=d->GetDCA();   
352         Double_t sigvert=d->GetSigmaVert();   
353         Double_t ptmax=0;
354         for(Int_t i=0;i<3;i++){
355           if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
356         }
357         fCosPHistLS[iPtBin]->Fill(cosp);
358         fDLenHistLS[iPtBin]->Fill(dlen);
359         fSumd02HistLS[iPtBin]->Fill(sumD02);
360         fSigVertHistLS[iPtBin]->Fill(sigvert);
361         fPtMaxHistLS[iPtBin]->Fill(ptmax);
362         fDCAHistLS[iPtBin]->Fill(dca);
363       }
364     }
365     if(unsetvtx) d->UnsetOwnPrimaryVtx();
366   }
367   //wei:normalized to the number of combinations (production)
368   //wei2:normalized to the number of  LS/OS (production)
369   //wei3:normalized to the number of  LS/OS (analysis)
370   //wei4:normalized to the number of combinations (analysis)
371   Float_t wei2=0;
372   if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
373   Float_t wei3=0;
374   if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)(nDplusLS+nDminusLS);
375   Float_t weiplus=1.,weiminus=1.;
376   Float_t wei4plus=1.,wei4minus=1.;
377   //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
378   if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
379   if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
380   if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
381   if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);
382
383   fMassHistLS[index]->Add(histLSPlus,weiplus);
384   fMassHistLS[index]->Add(histLSMinus,weiminus);
385   fMassHistLS[index+2]->Add(histLSPlus,wei2);
386   fMassHistLS[index+2]->Add(histLSMinus,wei2);
387   fMassHistLS[index+3]->Add(histLSPlus,wei3);
388   fMassHistLS[index+3]->Add(histLSMinus,wei3);
389   fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
390   fMassHistLS[index+4]->Add(histLSMinus,wei4minus);
391
392   delete histLSPlus;histLSPlus=0;
393   delete histLSMinus;histLSMinus=0;
394   
395   if(fDebug>1) printf("LS analysis terminated\n");  
396 }
397
398 //__________________________________________
399 void AliAnalysisTaskSEDplus::Init(){
400   //
401   // Initialization
402   //
403   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
404   
405   //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
406   fListCuts=new TList();
407   fListCuts->SetOwner();
408   AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
409   analysis->SetName("AnalysisCuts");
410   
411   fListCuts->Add(analysis);
412   PostData(2,fListCuts);
413   
414   return;
415 }
416
417 //________________________________________________________________________
418 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
419 {
420   // Create the output container
421   //
422   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
423
424   // Several histograms are more conveniently managed in a TList
425   fOutput = new TList();
426   fOutput->SetOwner();
427   fOutput->SetName("OutputHistos");
428
429   TString hisname;
430   Int_t index=0;
431   Int_t nbins=GetNBinsHistos();
432   fHistCentrality[0]=new TH2F("hCentrMult","centrality",100,0.5,30000.5,40,0.,100.);
433   fHistCentrality[1]=new TH2F("hCentrMult(selectedCent)","centrality(selectedCent)",100,0.5,30000.5,40,0.,100.);
434   fHistCentrality[2]=new TH2F("hCentrMult(OutofCent)","centrality(OutofCent)",100,0.5,30000.5,40,0.,100.);
435   for(Int_t i=0;i<3;i++){
436     fHistCentrality[i]->Sumw2();
437     fOutput->Add(fHistCentrality[i]);
438   }
439   for(Int_t i=0;i<fNPtBins;i++){
440
441     index=GetHistoIndex(i);
442
443     hisname.Form("hMassPt%d",i);
444     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
445     fMassHist[index]->Sumw2();
446     hisname.Form("hCosPAllPt%d",i);
447     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
448     fCosPHist[index]->Sumw2();
449     hisname.Form("hDLenAllPt%d",i);
450     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
451     fDLenHist[index]->Sumw2();
452     hisname.Form("hSumd02AllPt%d",i);
453     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
454     fSumd02Hist[index]->Sumw2();
455     hisname.Form("hSigVertAllPt%d",i);
456     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
457     fSigVertHist[index]->Sumw2();
458     hisname.Form("hPtMaxAllPt%d",i);
459     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
460     fPtMaxHist[index]->Sumw2();
461     hisname.Form("hPtKPt%d",i);
462     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
463     fPtKHist[index]->Sumw2();
464     hisname.Form("hPtpi1Pt%d",i);
465     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
466     fPtpi1Hist[index]->Sumw2();
467     hisname.Form("hPtpi2Pt%d",i);
468     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
469     fPtpi2Hist[index]->Sumw2();
470     hisname.Form("hDCAAllPt%d",i);
471     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
472     fDCAHist[index]->Sumw2();
473
474     hisname.Form("hDLxyPt%d",i);
475     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
476     fDLxy[index]->Sumw2();
477     hisname.Form("hCosxyPt%d",i);
478     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
479     fCosxy[index]->Sumw2();
480     hisname.Form("hDLxyPt%dTC",i);
481     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
482     fDLxyTC[index]->Sumw2();
483     hisname.Form("hCosxyPt%dTC",i);
484     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
485     fCosxyTC[index]->Sumw2();
486    
487     hisname.Form("hMassPt%dTC",i);
488     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
489     fMassHistTC[index]->Sumw2();
490     hisname.Form("hMassPt%dTCPlus",i);
491     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
492     fMassHistTCPlus[index]->Sumw2();
493     hisname.Form("hMassPt%dTCMinus",i);
494     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
495     fMassHistTCMinus[index]->Sumw2();
496
497
498        
499     index=GetSignalHistoIndex(i);    
500     hisname.Form("hSigPt%d",i);
501     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
502     fMassHist[index]->Sumw2();
503     hisname.Form("hCosPSigPt%d",i);
504     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
505     fCosPHist[index]->Sumw2();
506     hisname.Form("hDLenSigPt%d",i);
507     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
508     fDLenHist[index]->Sumw2();
509     hisname.Form("hSumd02SigPt%d",i);
510     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
511     fSumd02Hist[index]->Sumw2();
512     hisname.Form("hSigVertSigPt%d",i);
513     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
514     fSigVertHist[index]->Sumw2();
515     hisname.Form("hPtMaxSigPt%d",i);
516     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
517     fPtMaxHist[index]->Sumw2();  
518     hisname.Form("hPtKSigPt%d",i);  
519     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
520     fPtKHist[index]->Sumw2();
521     hisname.Form("hPtpi1SigPt%d",i);
522     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
523     fPtpi1Hist[index]->Sumw2();
524     hisname.Form("hPtpi2SigPt%d",i);
525     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
526     fPtpi2Hist[index]->Sumw2();
527
528     hisname.Form("hDCASigPt%d",i);
529     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
530     fDCAHist[index]->Sumw2();    
531
532     hisname.Form("hDLxySigPt%d",i);
533     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
534     fDLxy[index]->Sumw2();
535     hisname.Form("hCosxySigPt%d",i);
536     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
537     fCosxy[index]->Sumw2();
538     hisname.Form("hDLxySigPt%dTC",i);
539     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
540     fDLxyTC[index]->Sumw2();
541     hisname.Form("hCosxySigPt%dTC",i);
542     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
543     fCosxyTC[index]->Sumw2();
544     hisname.Form("hSigPt%dTC",i);
545     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
546     fMassHistTC[index]->Sumw2();
547     hisname.Form("hSigPt%dTCPlus",i);
548     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
549     fMassHistTCPlus[index]->Sumw2();
550     hisname.Form("hSigPt%dTCMinus",i);
551     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
552     fMassHistTCMinus[index]->Sumw2();
553
554
555     index=GetBackgroundHistoIndex(i); 
556     hisname.Form("hBkgPt%d",i);
557     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
558     fMassHist[index]->Sumw2();
559     hisname.Form("hCosPBkgPt%d",i);
560     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
561     fCosPHist[index]->Sumw2();
562     hisname.Form("hDLenBkgPt%d",i);
563     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
564     fDLenHist[index]->Sumw2();
565     hisname.Form("hSumd02BkgPt%d",i);
566     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
567     fSumd02Hist[index]->Sumw2();
568     hisname.Form("hSigVertBkgPt%d",i);
569     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
570     fSigVertHist[index]->Sumw2();
571     hisname.Form("hPtMaxBkgPt%d",i);
572     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
573     fPtMaxHist[index]->Sumw2();
574     hisname.Form("hPtKBkgPt%d",i);  
575     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
576     fPtKHist[index]->Sumw2();
577     hisname.Form("hPtpi1BkgPt%d",i);
578     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
579     fPtpi1Hist[index]->Sumw2();
580     hisname.Form("hPtpi2BkgPt%d",i);
581     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
582     fPtpi2Hist[index]->Sumw2();
583     hisname.Form("hDCABkgPt%d",i);
584     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
585     fDCAHist[index]->Sumw2();
586
587     hisname.Form("hDLxyBkgPt%d",i);
588     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
589     fDLxy[index]->Sumw2();
590     hisname.Form("hCosxyBkgPt%d",i);
591     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
592     fCosxy[index]->Sumw2();
593     hisname.Form("hDLxyBkgPt%dTC",i);
594     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
595     fDLxyTC[index]->Sumw2();
596     hisname.Form("hCosxyBkgPt%dTC",i);
597     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
598     fCosxyTC[index]->Sumw2();
599   
600
601     hisname.Form("hBkgPt%dTC",i);
602     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
603     fMassHistTC[index]->Sumw2();
604     hisname.Form("hBkgPt%dTCPlus",i);
605     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
606     fMassHistTCPlus[index]->Sumw2();
607     hisname.Form("hBkgPt%dTCMinus",i);
608     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
609     fMassHistTCMinus[index]->Sumw2();
610   }
611     
612
613   for(Int_t i=0; i<3*fNPtBins; i++){
614     fOutput->Add(fMassHist[i]);
615     if(fCutsDistr){  
616       fOutput->Add(fCosPHist[i]);
617       fOutput->Add(fDLenHist[i]);
618       fOutput->Add(fSumd02Hist[i]);
619       fOutput->Add(fSigVertHist[i]);
620       fOutput->Add(fPtMaxHist[i]);
621       fOutput->Add(fPtKHist[i]);
622       fOutput->Add(fPtpi1Hist[i]);
623       fOutput->Add(fPtpi2Hist[i]);
624       fOutput->Add(fDCAHist[i]);
625       fOutput->Add(fDLxy[i]);
626       fOutput->Add(fDLxyTC[i]);
627       fOutput->Add(fCosxy[i]);
628       fOutput->Add(fCosxyTC[i]);
629     }  
630     fOutput->Add(fMassHistTC[i]);
631     fOutput->Add(fMassHistTCPlus[i]);
632     fOutput->Add(fMassHistTCMinus[i]);
633
634   }
635   
636   if(fCutsDistr){  
637     fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
638     fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
639     fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
640     for(Int_t i=0; i<3; i++){
641       fCorreld0Kd0pi[i]->Sumw2();
642       fOutput->Add(fCorreld0Kd0pi[i]);
643     }
644   }
645   
646   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
647   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
648   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents accepted");
649   fHistNEvents->GetXaxis()->SetBinLabel(3,"Rejected due to trigger");
650   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected pileup events");
651   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to centrality"); 
652   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vtxz");
653   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to Physics Sel");
654   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
655   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
656   fHistNEvents->GetXaxis()->SetBinLabel(10,"D+ after loose cuts");
657   fHistNEvents->GetXaxis()->SetBinLabel(11,"D+ after tight cuts");
658  
659   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);  
660   fHistNEvents->Sumw2();
661   fHistNEvents->SetMinimum(0);
662   fOutput->Add(fHistNEvents);
663
664   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
665   fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);  
666   fYVsPt=new TH3F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
667   fYVsPtTC=new TH3F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
668   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
669   fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
670   fPhiEtaCand=new TH2F("hPhiEtaCand","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
671   fPhiEtaCandSigReg=new TH2F("hPhiEtaCandSigReg","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
672   fSPDMult = new TH1F("hSPDMult", "Tracklets multiplicity; Tracklets ; Entries",200,0.,200.); 
673   fOutput->Add(fPtVsMass);
674   fOutput->Add(fPtVsMassTC);
675   fOutput->Add(fYVsPt);
676   fOutput->Add(fYVsPtTC);
677   fOutput->Add(fYVsPtSig);
678   fOutput->Add(fYVsPtSigTC);
679   fOutput->Add(fPhiEtaCand);
680   fOutput->Add(fPhiEtaCandSigReg);
681   fOutput->Add(fSPDMult);
682
683
684   //Counter for Normalization
685   TString normName="NormalizationCounter";
686   AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
687   if(cont)normName=(TString)cont->GetName();
688   fCounter = new AliNormalizationCounter(normName.Data());
689   fCounter->Init();
690
691   if(fDoLS) CreateLikeSignHistos();
692   if(fDoImpPar) CreateImpactParameterHistos();
693
694   if(fFillNtuple){
695     OpenFile(4); // 4 is the slot number of the ntuple
696    
697     
698     fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:Pt:charge:piddau0:piddau1:piddau2:Ptpi:PtK:Ptpi2:mompi:momK:mompi2:cosp:cospxy:DecLen:NormDecLen:DecLenXY:NormDecLenXY:InvMass:sigvert:d0Pi:d0K:d0Pi2:maxdca:ntracks:centr:RunNumber:BadDau");
699
700   }
701   
702   return;
703 }
704
705 //________________________________________________________________________
706 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
707 {
708   // Execute analysis for current event:
709   // heavy flavor candidates association to MC truth
710
711   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
712   
713   TClonesArray *array3Prong = 0;
714   TClonesArray *arrayLikeSign =0;
715   if(!aod && AODEvent() && IsStandardAOD()) {
716     // In case there is an AOD handler writing a standard AOD, use the AOD 
717     // event in memory rather than the input (ESD) event.    
718     aod = dynamic_cast<AliAODEvent*> (AODEvent());
719     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
720     // have to taken from the AOD event hold by the AliAODExtension
721     AliAODHandler* aodHandler = (AliAODHandler*) 
722       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
723     if(aodHandler->GetExtensions()) {
724       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
725       AliAODEvent *aodFromExt = ext->GetAOD();
726       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
727       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
728     }
729   } else if(aod) {
730     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
731     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
732   }
733
734   if(!aod || !array3Prong) {
735     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
736     return;
737   }
738   if(!arrayLikeSign && fDoLS) {
739     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
740     return;
741   }
742
743   // fix for temporary bug in ESDfilter 
744   // the AODs with null vertex pointer didn't pass the PhysSel
745   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
746
747   //Store the event in AliNormalizationCounter->To disable for Pb-Pb? Add a flag to disable it?
748   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
749
750   fHistNEvents->Fill(0); // count event
751   Int_t runNumber=aod->GetRunNumber();
752
753   //Event selection
754   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
755   Float_t ntracks=aod->GetNTracks();
756   Float_t evCentr=0;
757   if(fRDCutsAnalysis->GetUseCentrality()>0) evCentr=fRDCutsAnalysis->GetCentrality(aod);
758   fHistCentrality[0]->Fill(ntracks,evCentr);
759   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(2);
760   if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(3); 
761   if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(4);fHistCentrality[2]->Fill(ntracks,evCentr);}
762   if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
763   if(fRDCutsAnalysis->GetWhyRejection()==7)fHistNEvents->Fill(6);
764
765   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
766   //TString trigclass=aod->GetFiredTriggerClasses();
767   // Post the data already here  
768   PostData(1,fOutput);
769   if(!isEvSel)return;
770   Int_t tracklets=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
771   // printf("ntracklet===%d\n",tracklets);
772   fSPDMult->Fill(tracklets);
773
774   fHistCentrality[1]->Fill(ntracks,evCentr);
775   fHistNEvents->Fill(1);
776
777   TClonesArray *arrayMC=0;
778   AliAODMCHeader *mcHeader=0;
779
780   // AOD primary vertex
781   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
782   //    vtx1->Print();
783   //   TString primTitle = vtx1->GetTitle();
784   //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
785  
786   // load MC particles
787   if(fReadMC){
788     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
789     if(!arrayMC) {
790       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
791       return;
792     }
793     
794     // load MC header
795     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
796     if(!mcHeader) {
797       printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
798       return;
799     }
800   }
801   
802   Int_t n3Prong = array3Prong->GetEntriesFast();
803   //  printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
804   
805   Int_t nOS=0;
806   Int_t index;
807   Int_t pdgDgDplustoKpipi[3]={321,211,211};
808
809   if(fDoLS>1){//Normalizations for LS
810     for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
811       AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
812       if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
813         if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod))nOS++;
814       }
815     }
816   }else{//Standard analysis
817     // Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};//TO REMOVE
818     //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
819     Int_t nSelectedloose=0,nSelectedtight=0;
820     for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
821       AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
822       fHistNEvents->Fill(7);
823       if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
824         fHistNEvents->Fill(8);
825         continue;
826       }
827     
828       Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
829
830       if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
831
832       Double_t etaD=d->Eta();
833       Double_t phiD=d->Phi();
834       if(fEtaSelection!=0){
835         if(fEtaSelection==1 && etaD<0) continue;
836         if(fEtaSelection==-1 && etaD>0) continue;
837       }
838
839       Bool_t unsetvtx=kFALSE;
840       if(!d->GetOwnPrimaryVtx()){
841         d->SetOwnPrimaryVtx(vtx1);
842         unsetvtx=kTRUE;
843       }
844
845       Double_t ptCand = d->Pt();
846       Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);     
847     
848       Bool_t recVtx=kFALSE;
849       AliAODVertex *origownvtx=0x0;
850       if(fRDCutsAnalysis->GetIsPrimaryWithoutDaughters()){
851         if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());  
852         if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
853         else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
854       }
855     
856       Int_t labDp=-1;
857       Bool_t isPrimary=kTRUE;
858       Float_t pdgCode=-2;
859       Float_t trueImpParXY=0.;
860       if(fReadMC){
861         labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
862         if(labDp>=0){
863           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
864           if(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE;
865           pdgCode=TMath::Abs(partDp->GetPdgCode());
866           if(!isPrimary){
867             trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
868           }
869         }else{
870           pdgCode=-1;
871         }
872       }
873
874       Double_t invMass=d->InvMassDplus();
875       Double_t rapid=d->YDplus();
876       fYVsPt->Fill(ptCand,rapid,invMass);
877       if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid,invMass);nOS++;}
878       Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
879       if(isFidAcc){
880         fPtVsMass->Fill(invMass,ptCand);
881         if(passTightCuts){
882           fPtVsMassTC->Fill(invMass,ptCand);
883           fPhiEtaCand->Fill(etaD,phiD);
884           if(TMath::Abs(invMass-1.8696)<0.05) fPhiEtaCandSigReg->Fill(etaD,phiD);
885         }
886       }
887
888
889       Double_t  dlen=0,cosp=0,maxdca=0,sigvert=0,sumD02=0,ptmax=0,dlxy=0,cxy=0;
890       if(fCutsDistr||fFillNtuple||fDoImpPar){
891         dlen=d->DecayLength();
892         cosp=d->CosPointingAngle();
893         sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
894         maxdca=-9999.;
895         for(Int_t idau=0;idau<3;idau++) if(d->GetDCA(idau)>maxdca) maxdca=d->GetDCA(idau);
896         sigvert=d->GetSigmaVert();         
897         ptmax=0;
898         for(Int_t i=0;i<3;i++){
899           if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
900         }
901         dlxy=d->NormalizedDecayLengthXY();
902         cxy=d->CosPointingAngleXY();
903       }
904       Double_t impparXY=d->ImpParXY()*10000.;
905       Double_t arrayForSparse[6]={invMass,ptCand,impparXY,cosp,dlen,static_cast<Double_t>(tracklets)};
906       Double_t arrayForSparseTrue[6]={invMass,ptCand,trueImpParXY,cosp,dlen,static_cast<Double_t>(tracklets)};
907
908       //Ntuple
909       Float_t tmp[31];
910       if(fFillNtuple){
911         tmp[0]=pdgCode;
912         if(!isPrimary) tmp[0]+=5000.;
913         tmp[1]=d->Px();
914         tmp[2]=d->Py();
915         tmp[3]=d->Pz();
916         tmp[4]=d->Pt();
917         tmp[5]=d->GetCharge();
918         //      tmp[5]=fRDCutsAnalysis->GetPIDBitMask(d);         
919         tmp[6]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(0));
920         tmp[7]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(1));
921         tmp[8]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(2));
922         tmp[9]=d->PtProng(0);
923         tmp[10]=d->PtProng(1);
924         tmp[11]=d->PtProng(2);
925         tmp[12]=d->PProng(0);
926         tmp[13]=d->PProng(1);
927         tmp[14]=d->PProng(2);
928         tmp[15]=cosp;
929         tmp[16]=cxy;
930         tmp[17]=dlen;
931         tmp[18]=d->NormalizedDecayLength();
932         tmp[19]=d->DecayLengthXY();
933         tmp[20]=dlxy;
934         tmp[21]=d->InvMassDplus();
935         tmp[22]=sigvert;
936         tmp[23]=d->Getd0Prong(0);
937         tmp[24]=d->Getd0Prong(1);
938         tmp[25]=d->Getd0Prong(2);
939         tmp[26]=maxdca;
940         tmp[27]=ntracks;
941         tmp[28]=fRDCutsAnalysis->GetCentrality(aod);
942         tmp[29]=runNumber;
943         tmp[30]=d->HasBadDaughters();
944         fNtupleDplus->Fill(tmp);
945         PostData(4,fNtupleDplus);
946       }
947
948       //Fill histos
949       index=GetHistoIndex(iPtBin);
950       if(isFidAcc){
951         fHistNEvents->Fill(9);
952         nSelectedloose++;
953         fMassHist[index]->Fill(invMass);
954         if(fCutsDistr){   
955           fCosPHist[index]->Fill(cosp);
956           fDLenHist[index]->Fill(dlen);
957           fSumd02Hist[index]->Fill(sumD02);
958           fSigVertHist[index]->Fill(sigvert);
959           fPtMaxHist[index]->Fill(ptmax);
960           fPtKHist[index]->Fill(d->PtProng(1));
961           fPtpi1Hist[index]->Fill(d->PtProng(0));
962           fPtpi2Hist[index]->Fill(d->PtProng(2));
963           fDCAHist[index]->Fill(maxdca);
964           fDLxy[index]->Fill(dlxy);
965           fCosxy[index]->Fill(cxy);
966           fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
967                                   d->Getd0Prong(2)*d->Getd0Prong(1));
968         }
969         if(passTightCuts){ 
970           fHistNEvents->Fill(10);
971           nSelectedtight++;
972           fMassHistTC[index]->Fill(invMass);
973           if(fCutsDistr){  
974             fDLxyTC[index]->Fill(dlxy);
975             fCosxyTC[index]->Fill(cxy);
976           }
977           if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
978           else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
979           if(fDoImpPar){
980             fHistMassPtImpParTC[0]->Fill(arrayForSparse);
981           }
982         }
983       }
984       
985       if(fReadMC){
986         if(isFidAcc){
987           Int_t correlIndex=0;
988           if(labDp>=0) {
989             index=GetSignalHistoIndex(iPtBin);
990             correlIndex=1;
991             if(passTightCuts&&fDoImpPar){
992               if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
993               else{
994                 fHistMassPtImpParTC[2]->Fill(arrayForSparse);
995                 fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
996               }
997             }
998           }else{
999             index=GetBackgroundHistoIndex(iPtBin);
1000             correlIndex=2;
1001             if(passTightCuts&&fDoImpPar)fHistMassPtImpParTC[4]->Fill(arrayForSparse);
1002           }
1003           
1004           fMassHist[index]->Fill(invMass);
1005           if(fCutsDistr){
1006             Float_t fact=1.;
1007             Float_t factor[3]={1.,1.,1.};
1008             if(fUseStrangeness) fact=GetStrangenessWeights(d,arrayMC,factor);
1009             fCosPHist[index]->Fill(cosp,fact);
1010             fDLenHist[index]->Fill(dlen,fact);
1011             fDLxy[index]->Fill(dlxy);
1012             fCosxy[index]->Fill(cxy);
1013             
1014             Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2];
1015             fSumd02Hist[index]->Fill(sumd02s);
1016             fSigVertHist[index]->Fill(sigvert,fact);
1017             fPtMaxHist[index]->Fill(ptmax,fact);
1018             fPtKHist[index]->Fill(d->PtProng(1),fact);
1019             fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1020             fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1021             fDCAHist[index]->Fill(maxdca,fact);
1022             fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1023                                     d->Getd0Prong(2)*d->Getd0Prong(1));
1024           }
1025           if(passTightCuts){
1026             fMassHistTC[index]->Fill(invMass);        
1027             if(fCutsDistr){
1028               fDLxyTC[index]->Fill(dlxy);
1029               fCosxyTC[index]->Fill(cxy);
1030             }         
1031             if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1032             else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1033           }
1034         }else{//outside fidAcc
1035           if(labDp>=0){
1036             fYVsPtSig->Fill(ptCand,rapid, invMass);
1037             if(passTightCuts)fYVsPtSigTC->Fill(ptCand,rapid, invMass);
1038           }
1039         }
1040       }//readmc
1041     
1042       if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
1043       
1044       if(unsetvtx) d->UnsetOwnPrimaryVtx();
1045     }
1046     fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1047     fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1048   }
1049   //start LS analysis
1050   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1051   
1052   PostData(1,fOutput); 
1053   PostData(2,fListCuts);
1054   PostData(3,fCounter);    
1055   return;
1056 }
1057
1058 //________________________________________________________________________
1059 void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){
1060   // Histos for Like Sign bckground
1061
1062   TString hisname;
1063   Int_t indexLS=0;
1064   Int_t index=0;
1065   Int_t nbins=GetNBinsHistos();
1066   for(Int_t i=0;i<fNPtBins;i++){
1067
1068     index=GetHistoIndex(i);
1069     indexLS=GetLSHistoIndex(i);
1070
1071     hisname.Form("hLSPt%dLC",i);
1072     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1073     fMassHistLS[indexLS]->Sumw2();
1074     hisname.Form("hLSPt%dTC",i);
1075     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1076     fMassHistLSTC[indexLS]->Sumw2();
1077
1078     hisname.Form("hCosPAllPt%dLS",i);
1079     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1080     fCosPHistLS[index]->Sumw2();
1081     hisname.Form("hDLenAllPt%dLS",i);
1082     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1083     fDLenHistLS[index]->Sumw2();
1084     hisname.Form("hSumd02AllPt%dLS",i);
1085     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1086     fSumd02HistLS[index]->Sumw2();
1087     hisname.Form("hSigVertAllPt%dLS",i);
1088     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1089     fSigVertHistLS[index]->Sumw2();
1090     hisname.Form("hPtMaxAllPt%dLS",i);
1091     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1092     fPtMaxHistLS[index]->Sumw2();
1093     hisname.Form("hDCAAllPt%dLS",i);
1094     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1095     fDCAHistLS[index]->Sumw2();    
1096
1097     index=GetSignalHistoIndex(i);    
1098     indexLS++;
1099  
1100     hisname.Form("hLSPt%dLCnw",i);
1101     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1102     fMassHistLS[indexLS]->Sumw2();
1103     hisname.Form("hLSPt%dTCnw",i);
1104     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1105     fMassHistLSTC[indexLS]->Sumw2();
1106
1107     hisname.Form("hCosPSigPt%dLS",i);
1108     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1109     fCosPHistLS[index]->Sumw2();
1110     hisname.Form("hDLenSigPt%dLS",i);
1111     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1112     fDLenHistLS[index]->Sumw2();
1113     hisname.Form("hSumd02SigPt%dLS",i);
1114     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1115     fSumd02HistLS[index]->Sumw2();
1116     hisname.Form("hSigVertSigPt%dLS",i);
1117     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1118     fSigVertHistLS[index]->Sumw2();
1119     hisname.Form("hPtMaxSigPt%dLS",i);
1120     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1121     fPtMaxHistLS[index]->Sumw2();
1122     hisname.Form("hDCASigPt%dLS",i);
1123     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1124     fDCAHistLS[index]->Sumw2();
1125
1126     index=GetBackgroundHistoIndex(i); 
1127     indexLS++;
1128
1129     hisname.Form("hLSPt%dLCntrip",i);
1130     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1131     fMassHistLS[indexLS]->Sumw2();
1132     hisname.Form("hLSPt%dTCntrip",i);
1133     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1134     fMassHistLSTC[indexLS]->Sumw2();
1135
1136     hisname.Form("hCosPBkgPt%dLS",i);
1137     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1138     fCosPHistLS[index]->Sumw2();
1139     hisname.Form("hDLenBkgPt%dLS",i);
1140     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1141     fDLenHistLS[index]->Sumw2();
1142     hisname.Form("hSumd02BkgPt%dLS",i);
1143     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1144     fSumd02HistLS[index]->Sumw2();
1145     hisname.Form("hSigVertBkgPt%dLS",i);
1146     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1147     fSigVertHistLS[index]->Sumw2();
1148     hisname.Form("hPtMaxBkgPt%dLS",i);
1149     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1150     fPtMaxHistLS[index]->Sumw2();
1151     hisname.Form("hDCABkgPt%dLS",i);
1152     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1153     fDCAHistLS[index]->Sumw2();
1154
1155     indexLS++;
1156     hisname.Form("hLSPt%dLCntripsinglecut",i);
1157     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1158     fMassHistLS[indexLS]->Sumw2();
1159     hisname.Form("hLSPt%dTCntripsinglecut",i);
1160     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1161     fMassHistLSTC[indexLS]->Sumw2();
1162
1163     indexLS++;
1164     hisname.Form("hLSPt%dLCspc",i);
1165     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1166     fMassHistLS[indexLS]->Sumw2();
1167     hisname.Form("hLSPt%dTCspc",i);
1168     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1169     fMassHistLSTC[indexLS]->Sumw2();
1170   }
1171
1172   for(Int_t i=0; i<3*fNPtBins; i++){
1173     fOutput->Add(fCosPHistLS[i]);
1174     fOutput->Add(fDLenHistLS[i]);
1175     fOutput->Add(fSumd02HistLS[i]);
1176     fOutput->Add(fSigVertHistLS[i]);
1177     fOutput->Add(fPtMaxHistLS[i]);  
1178     fOutput->Add(fDCAHistLS[i]);  
1179   }
1180   for(Int_t i=0; i<5*fNPtBins; i++){
1181     fOutput->Add(fMassHistLS[i]);
1182     fOutput->Add(fMassHistLSTC[i]);
1183   }
1184 }
1185
1186 //________________________________________________________________________
1187 void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){
1188   // Histos for impact paramter study
1189
1190   Int_t nmassbins=GetNBinsHistos();
1191   Double_t maxmult;
1192   if(fSystem==1) maxmult=5000;
1193   else maxmult=200;
1194   Int_t nbins[6]={nmassbins,200,fNImpParBins,5,50,100};
1195   Double_t xmin[6]={fLowmasslimit,0.,fLowerImpPar,0.95,0.,-0.5};
1196   Double_t xmax[6]={fUpmasslimit,40.,fHigherImpPar,1.,1.,maxmult};
1197
1198   fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1199                                         "Mass vs. pt vs.imppar - All",
1200                                         6,nbins,xmin,xmax);
1201   fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1202                                         "Mass vs. pt vs.imppar - promptD",
1203                                         6,nbins,xmin,xmax);
1204   fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1205                                         "Mass vs. pt vs.imppar - DfromB",
1206                                         6,nbins,xmin,xmax);
1207   fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1208                                         "Mass vs. pt vs.true imppar -DfromB",
1209                                         6,nbins,xmin,xmax);
1210   fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1211                                         "Mass vs. pt vs.imppar - backgr.",
1212                                         6,nbins,xmin,xmax);
1213   for(Int_t i=0; i<5;i++){
1214     fOutput->Add(fHistMassPtImpParTC[i]);
1215   }
1216 }
1217
1218 //________________________________________________________________________
1219 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1220 {
1221   // Terminate analysis
1222   //
1223   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1224
1225   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1226   if (!fOutput) {     
1227     printf("ERROR: fOutput not available\n");
1228     return;
1229   }
1230
1231   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1232   if(fHistNEvents){
1233     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1234   }else{
1235     printf("ERROR: fHistNEvents not available\n");
1236     return;
1237   }
1238
1239   return;
1240 }
1241 //_________________________________________________________________________________________________
1242 Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, const AliAODMCParticle *mcPartCandidate) const {               
1243   //
1244   // checking whether the mother of the particles come from a charm or a bottom quark
1245   //
1246         
1247   Int_t pdgGranma = 0;
1248   Int_t mother = 0;
1249   mother = mcPartCandidate->GetMother();
1250   Int_t istep = 0;
1251   Int_t abspdgGranma =0;
1252   Bool_t isFromB=kFALSE;
1253   Bool_t isQuarkFound=kFALSE;
1254   while (mother >0 ){
1255     istep++;
1256     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1257     if (mcGranma){
1258       pdgGranma = mcGranma->GetPdgCode();
1259       abspdgGranma = TMath::Abs(pdgGranma);
1260       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1261         isFromB=kTRUE;
1262       }
1263       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1264       mother = mcGranma->GetMother();
1265     }else{
1266       AliError("Failed casting the mother particle!");
1267       break;
1268     }
1269   }
1270   
1271   if(isFromB) return 5;
1272   else return 4;
1273 }
1274 //_________________________________________________________________________________________________
1275 Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1276   // true impact parameter calculation
1277
1278   Double_t vtxTrue[3];
1279   mcHeader->GetVertex(vtxTrue);
1280   Double_t origD[3];
1281   partDp->XvYvZv(origD);
1282   Short_t charge=partDp->Charge();
1283   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1284   for(Int_t iDau=0; iDau<3; iDau++){
1285     pXdauTrue[iDau]=0.;
1286     pYdauTrue[iDau]=0.;
1287     pZdauTrue[iDau]=0.;
1288   }
1289
1290   Int_t nDau=partDp->GetNDaughters();
1291   Int_t labelFirstDau = partDp->GetDaughter(0); 
1292   if(nDau==3){
1293     for(Int_t iDau=0; iDau<3; iDau++){
1294       Int_t ind = labelFirstDau+iDau;
1295       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1296       if(!part){
1297         AliError("Daughter particle not found in MC array");
1298         return 99999.;
1299       } 
1300       pXdauTrue[iDau]=part->Px();
1301       pYdauTrue[iDau]=part->Py();
1302       pZdauTrue[iDau]=part->Pz();
1303     }
1304   }else if(nDau==2){
1305     Int_t theDau=0;
1306     for(Int_t iDau=0; iDau<2; iDau++){
1307       Int_t ind = labelFirstDau+iDau;
1308       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1309       if(!part){
1310         AliError("Daughter particle not found in MC array");
1311         return 99999.;
1312       } 
1313       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1314       if(pdgCode==211 || pdgCode==321){
1315         pXdauTrue[theDau]=part->Px();
1316         pYdauTrue[theDau]=part->Py();
1317         pZdauTrue[theDau]=part->Pz();
1318         ++theDau;
1319       }else{
1320         Int_t nDauRes=part->GetNDaughters();
1321         if(nDauRes==2){
1322           Int_t labelFirstDauRes = part->GetDaughter(0);        
1323           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1324             Int_t indDR = labelFirstDauRes+iDauRes;
1325             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1326             if(!partDR){
1327               AliError("Daughter particle not found in MC array");
1328               return 99999.;
1329             } 
1330             
1331             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1332             if(pdgCodeDR==211 || pdgCodeDR==321){
1333               pXdauTrue[theDau]=partDR->Px();
1334               pYdauTrue[theDau]=partDR->Py();
1335               pZdauTrue[theDau]=partDR->Pz();
1336               ++theDau;
1337             }
1338           }
1339         }
1340       }
1341     }
1342     if(theDau!=3){
1343       AliError("Wrong number of decay prongs");
1344       return 99999.;
1345     }
1346   }
1347
1348   Double_t d0dummy[3]={0.,0.,0.};
1349   AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1350   return aodDplusMC.ImpParXY();
1351
1352 }
1353 //_________________________________________________________________________________________________
1354 Float_t AliAnalysisTaskSEDplus::GetStrangenessWeights(const AliAODRecoDecayHF3Prong* d, TClonesArray* arrayMC, Float_t factor[3]) const {
1355   // Computes weights to adapt strangeness in MC to data
1356
1357   for(Int_t iprong=0;iprong<3;iprong++){
1358     factor[iprong]=1;
1359     AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1360     Int_t labd= trad->GetLabel();
1361     if(labd>=0){
1362       AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1363       if(dau){
1364         Int_t labm = dau->GetMother();
1365         if(labm>=0){
1366           AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1367           if(mot){
1368             if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1369               if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1370               else factor[iprong]=1./.6;
1371               //          fNentries->Fill(12);
1372             }
1373             if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1374               factor[iprong]=1./0.25;
1375               //                  fNentries->Fill(13);
1376             }//if 3122
1377           }//if(mot)
1378         }//if labm>0
1379       }//if(dau)
1380     }//if labd>=0
1381   }//prong loop
1382
1383   Float_t fact=1.;
1384   for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1385   return fact;
1386
1387 }