Updates for proper treatment of fQuality flag in centrality selection and for usage...
[u/mrichter/AliRoot.git] / PWG3 / 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 <TH1F.h>
34 #include <TDatabasePDG.h>
35
36 #include "AliAnalysisManager.h"
37 #include "AliRDHFCutsDplustoKpipi.h"
38 #include "AliAODHandler.h"
39 #include "AliAODEvent.h"
40 #include "AliAODVertex.h"
41 #include "AliAODTrack.h"
42 #include "AliAODMCHeader.h"
43 #include "AliAODMCParticle.h"
44 #include "AliAODRecoDecayHF3Prong.h"
45 #include "AliAnalysisVertexingHF.h"
46 #include "AliAnalysisTaskSE.h"
47 #include "AliAnalysisTaskSEDplus.h"
48 #include "AliNormalizationCounter.h"
49
50 ClassImp(AliAnalysisTaskSEDplus)
51
52
53 //________________________________________________________________________
54 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
55 AliAnalysisTaskSE(),
56   fOutput(0), 
57   fHistNEvents(0),
58   fPtVsMass(0),
59   fPtVsMassTC(0),
60   fYVsPt(0),
61   fYVsPtTC(0),
62   fYVsPtSig(0),
63   fYVsPtSigTC(0),
64   fNtupleDplus(0),
65   fUpmasslimit(1.965),
66   fLowmasslimit(1.765),
67   fNPtBins(0),
68   fBinWidth(0.002),
69   fListCuts(0),
70   fRDCutsProduction(0),
71   fRDCutsAnalysis(0),
72   fCounter(0),
73   fFillNtuple(kFALSE),
74   fReadMC(kFALSE),
75   fUseStrangeness(kFALSE),
76   fUseBit(kTRUE),
77   fDoLS(kFALSE)
78 {
79    // Default constructor
80 }
81
82 //________________________________________________________________________
83 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple):
84 AliAnalysisTaskSE(name),
85 fOutput(0),
86 fHistNEvents(0),
87 fPtVsMass(0),
88 fPtVsMassTC(0),
89 fYVsPt(0),
90 fYVsPtTC(0),
91 fYVsPtSig(0),
92 fYVsPtSigTC(0),
93 fNtupleDplus(0),
94 fUpmasslimit(1.965),
95 fLowmasslimit(1.765),
96 fNPtBins(0),
97 fBinWidth(0.002),
98 fListCuts(0),
99 fRDCutsProduction(dpluscutsprod),
100 fRDCutsAnalysis(dpluscutsana),
101 fCounter(0),
102 fFillNtuple(fillNtuple),
103 fReadMC(kFALSE),
104 fUseStrangeness(kFALSE),
105 fUseBit(kTRUE),
106 fDoLS(kFALSE)
107 {
108   // 
109   // Standrd constructor
110   //
111   //Double_t ptlim[5]={0.,2.,3.,5,9999999.};
112    //SetPtBinLimit(5, ptlim);
113   SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
114   // Default constructor
115    // Output slot #1 writes into a TList container
116   DefineOutput(1,TList::Class());  //My private output
117  // Output slot #2 writes cut to private output
118   //  DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
119   DefineOutput(2,TList::Class());
120 // Output slot #3 writes cut to private output
121   DefineOutput(3,AliNormalizationCounter::Class());
122
123   if(fFillNtuple){
124     // Output slot #4 writes into a TNtuple container
125     DefineOutput(4,TNtuple::Class());  //My private output
126   }
127 }
128
129 //________________________________________________________________________
130 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
131 {
132   //
133   // Destructor
134   //
135   if (fOutput) {
136     delete fOutput;
137     fOutput = 0;
138   }
139   if(fHistNEvents){
140     delete fHistNEvents;
141     fHistNEvents=0;
142   }  
143   for(Int_t i=0;i<3;i++){
144     if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;}
145   }
146   for(Int_t i=0;i<3*fNPtBins;i++){
147     if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
148     if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
149     if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
150     if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
151     if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
152     if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
153     if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
154     if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
155     if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
156     if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
157
158     if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
159     if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
160     if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
161     if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
162     if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
163     if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
164     if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
165     if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
166   }
167   if(fPtVsMass){
168     delete fPtVsMass;
169     fPtVsMass=0;
170   }
171   if(fPtVsMassTC){
172     delete fPtVsMassTC;
173     fPtVsMassTC=0;
174   }
175   if(fYVsPt){
176     delete fYVsPt;
177     fYVsPt=0;
178   }
179   if(fYVsPtTC){
180     delete fYVsPtTC;
181     fYVsPtTC=0;
182   }
183   if(fYVsPtSig){
184     delete fYVsPtSig;
185     fYVsPtSig=0;
186   }
187   if(fYVsPtSigTC){
188     delete fYVsPtSigTC;
189     fYVsPtSigTC=0;
190   }
191   if(fNtupleDplus){
192     delete fNtupleDplus;
193     fNtupleDplus=0;
194   }
195   if (fListCuts) {
196     delete fListCuts;
197     fListCuts = 0;
198   }
199   if(fRDCutsProduction){
200     delete fRDCutsProduction;
201     fRDCutsProduction = 0;
202   }
203   if(fRDCutsAnalysis){
204     delete fRDCutsAnalysis;
205     fRDCutsAnalysis = 0;
206   }
207
208   if(fCounter){
209     delete fCounter;
210     fCounter = 0;
211   }
212
213
214 }  
215 //_________________________________________________________________
216 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
217   // set invariant mass limits
218   Float_t bw=GetBinWidth();
219   fUpmasslimit = 1.865+range;
220   fLowmasslimit = 1.865-range;
221   SetBinWidth(bw);
222 }
223 //_________________________________________________________________
224 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
225   // set invariant mass limits
226   if(uplimit>lowlimit)
227     {
228       Float_t bw=GetBinWidth();
229       fUpmasslimit = lowlimit;
230       fLowmasslimit = uplimit;
231       SetBinWidth(bw);
232     }
233 }
234 //________________________________________________________________________
235 void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Float_t* lim){
236   // define pt bins for analysis
237   if(n>kMaxPtBins){
238     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
239     fNPtBins=kMaxPtBins;
240     fArrayBinLimits[0]=0.;
241     fArrayBinLimits[1]=2.;
242     fArrayBinLimits[2]=3.;
243     fArrayBinLimits[3]=5.;
244     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
245   }else{
246     fNPtBins=n-1;
247     fArrayBinLimits[0]=lim[0];
248     for(Int_t i=1; i<fNPtBins+1; i++) 
249       if(lim[i]>fArrayBinLimits[i-1]){
250         fArrayBinLimits[i]=lim[i];
251       }
252       else {
253         fArrayBinLimits[i]=fArrayBinLimits[i-1];
254       }
255     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
256   }
257   if(fDebug > 1){
258     printf("Number of Pt bins = %d\n",fNPtBins);
259     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);    
260   }
261 }
262 //________________________________________________________________
263 void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
264   Float_t width=w;
265   Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
266   Int_t missingbins=4-nbins%4;
267   nbins=nbins+missingbins;
268   width=(fUpmasslimit-fLowmasslimit)/nbins;
269   if(missingbins!=0){
270     printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
271   }
272   else{
273     if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
274   }
275   fBinWidth=width;
276 }
277 //_________________________________________________________________
278 Double_t  AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
279   // get pt bin limit
280   if(ibin>fNPtBins)return -1;
281   return fArrayBinLimits[ibin];
282
283 //_________________________________________________________________
284 Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
285   return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
286 }
287 //_________________________________________________________________
288 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
289   //
290   //
291   // Fill the Like Sign histograms
292   //
293   printf("started LS\n");
294   Bool_t fPlotVariables=kTRUE;
295   //histograms for like sign
296   Int_t nbins=GetNBinsHistos();;
297   TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
298   TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
299   
300   Int_t nPosTrks=0,nNegTrks=0;
301   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
302   Int_t nDplusLS=0;
303   Int_t nDminusLS=0;
304   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
305   Int_t index=0; 
306   
307   // loop over like sign candidates
308   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
309     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
310     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
311     Bool_t unsetvtx=kFALSE;
312     if(!d->GetOwnPrimaryVtx()) {
313       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
314       unsetvtx=kTRUE;
315     }
316     
317     //    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){
318     fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
319     Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
320     if(passTightCuts>0){
321       
322       Double_t ptCand = d->Pt();
323       Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
324       if(iPtBin<0)continue;
325       
326       Int_t sign= d->GetCharge();
327       Float_t invMass = d->InvMassDplus();
328       
329       index=GetLSHistoIndex(iPtBin);
330       fMassHistLS[index+1]->Fill(invMass);
331       if(sign>0){
332         histLSPlus->Fill(invMass);
333         nDplusLS++;
334       }else{
335         histLSMinus->Fill(invMass);
336         nDminusLS++;
337       }
338       if(fPlotVariables){
339         Double_t dlen=d->DecayLength();
340         Double_t cosp=d->CosPointingAngle();
341         Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
342         Double_t dca=d->GetDCA();   
343         Double_t sigvert=d->GetSigmaVert();   
344         Double_t ptmax=0;
345         for(Int_t i=0;i<3;i++){
346           if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
347         }
348         fCosPHistLS[iPtBin]->Fill(cosp);
349         fDLenHistLS[iPtBin]->Fill(dlen);
350         fSumd02HistLS[iPtBin]->Fill(sumD02);
351         fSigVertHistLS[iPtBin]->Fill(sigvert);
352         fPtMaxHistLS[iPtBin]->Fill(ptmax);
353         fDCAHistLS[iPtBin]->Fill(dca);
354       }
355     }
356     if(unsetvtx) d->UnsetOwnPrimaryVtx();
357   }
358   //wei:normalized to the number of combinations (production)
359   //wei2:normalized to the number of  LS/OS (production)
360   //wei3:normalized to the number of  LS/OS (analysis)
361   //wei4:normalized to the number of combinations (analysis)
362   Float_t wei2=0;
363   if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
364   Float_t wei3=0;
365   if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
366   Float_t weiplus=1.,weiminus=1.;
367   Float_t wei4plus=1.,wei4minus=1.;
368   //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
369   if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
370   if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
371   if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
372   if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);
373
374   fMassHistLS[index]->Add(histLSPlus,weiplus);
375   fMassHistLS[index]->Add(histLSMinus,weiminus);
376   fMassHistLS[index+2]->Add(histLSPlus,wei2);
377   fMassHistLS[index+2]->Add(histLSMinus,wei2);
378   fMassHistLS[index+3]->Add(histLSPlus,wei3);
379   fMassHistLS[index+3]->Add(histLSMinus,wei3);
380   fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
381   fMassHistLS[index+4]->Add(histLSMinus,wei4minus);
382
383   delete histLSPlus;histLSPlus=0;
384   delete histLSMinus;histLSMinus=0;
385   
386   if(fDebug>1) printf("LS analysis terminated\n");  
387 }
388
389 //__________________________________________
390 void AliAnalysisTaskSEDplus::Init(){
391   //
392   // Initialization
393   //
394   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
395   
396   //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
397   fListCuts=new TList();
398   AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction);
399   production->SetName("ProductionCuts");
400   AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
401   analysis->SetName("AnalysisCuts");
402   
403   fListCuts->Add(production);
404   fListCuts->Add(analysis);
405   PostData(2,fListCuts);
406   
407   return;
408 }
409
410 //________________________________________________________________________
411 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
412 {
413   // Create the output container
414   //
415   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
416
417   // Several histograms are more conveniently managed in a TList
418   fOutput = new TList();
419   fOutput->SetOwner();
420   fOutput->SetName("OutputHistos");
421
422   TString hisname;
423   Int_t index=0;
424   Int_t indexLS=0;
425   Int_t nbins=GetNBinsHistos();
426   fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5);
427   fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5);
428   fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5);
429   for(Int_t i=0;i<3;i++){
430     fHistCentrality[i]->Sumw2();
431     fOutput->Add(fHistCentrality[i]);
432   }
433   for(Int_t i=0;i<fNPtBins;i++){
434
435     index=GetHistoIndex(i);
436     indexLS=GetLSHistoIndex(i);
437
438     hisname.Form("hMassPt%d",i);
439     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
440     fMassHist[index]->Sumw2();
441     hisname.Form("hCosPAllPt%d",i);
442     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
443     fCosPHist[index]->Sumw2();
444     hisname.Form("hDLenAllPt%d",i);
445     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
446     fDLenHist[index]->Sumw2();
447     hisname.Form("hSumd02AllPt%d",i);
448     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
449     fSumd02Hist[index]->Sumw2();
450     hisname.Form("hSigVertAllPt%d",i);
451     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
452     fSigVertHist[index]->Sumw2();
453     hisname.Form("hPtMaxAllPt%d",i);
454     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
455     fPtMaxHist[index]->Sumw2();
456
457     hisname.Form("hDCAAllPt%d",i);
458     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
459     fDCAHist[index]->Sumw2();
460
461
462
463     hisname.Form("hMassPt%dTC",i);
464     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
465     fMassHistTC[index]->Sumw2();
466     hisname.Form("hMassPt%dTCPlus",i);
467     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
468     fMassHistTCPlus[index]->Sumw2();
469     hisname.Form("hMassPt%dTCMinus",i);
470     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
471     fMassHistTCMinus[index]->Sumw2();
472
473
474
475     
476     hisname.Form("hCosPAllPt%dLS",i);
477     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
478     fCosPHistLS[index]->Sumw2();
479     hisname.Form("hDLenAllPt%dLS",i);
480     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
481     fDLenHistLS[index]->Sumw2();
482     hisname.Form("hSumd02AllPt%dLS",i);
483     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
484     fSumd02HistLS[index]->Sumw2();
485     hisname.Form("hSigVertAllPt%dLS",i);
486     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
487     fSigVertHistLS[index]->Sumw2();
488     hisname.Form("hPtMaxAllPt%dLS",i);
489     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
490     fPtMaxHistLS[index]->Sumw2();
491     
492     hisname.Form("hDCAAllPt%dLS",i);
493     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
494     fDCAHistLS[index]->Sumw2();
495     
496     hisname.Form("hLSPt%dLC",i);
497     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
498     fMassHistLS[indexLS]->Sumw2();
499     
500     hisname.Form("hLSPt%dTC",i);
501     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
502     fMassHistLSTC[indexLS]->Sumw2();
503
504
505     
506     index=GetSignalHistoIndex(i);    
507     indexLS++;
508     hisname.Form("hSigPt%d",i);
509     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
510     fMassHist[index]->Sumw2();
511     hisname.Form("hCosPSigPt%d",i);
512     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
513     fCosPHist[index]->Sumw2();
514     hisname.Form("hDLenSigPt%d",i);
515     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
516     fDLenHist[index]->Sumw2();
517     hisname.Form("hSumd02SigPt%d",i);
518     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
519     fSumd02Hist[index]->Sumw2();
520     hisname.Form("hSigVertSigPt%d",i);
521     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
522     fSigVertHist[index]->Sumw2();
523     hisname.Form("hPtMaxSigPt%d",i);
524     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
525     fPtMaxHist[index]->Sumw2();    
526
527     hisname.Form("hDCASigPt%d",i);
528     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
529     fDCAHist[index]->Sumw2();    
530
531
532     hisname.Form("hSigPt%dTC",i);
533     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
534     fMassHistTC[index]->Sumw2();
535     hisname.Form("hSigPt%dTCPlus",i);
536     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
537     fMassHistTCPlus[index]->Sumw2();
538     hisname.Form("hSigPt%dTCMinus",i);
539     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
540     fMassHistTCMinus[index]->Sumw2();
541
542     hisname.Form("hLSPt%dLCnw",i);
543     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
544     fMassHistLS[indexLS]->Sumw2();
545     hisname.Form("hLSPt%dTCnw",i);
546     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
547     fMassHistLSTC[indexLS]->Sumw2();
548
549
550     
551     hisname.Form("hCosPSigPt%dLS",i);
552     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
553     fCosPHistLS[index]->Sumw2();
554     hisname.Form("hDLenSigPt%dLS",i);
555     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
556     fDLenHistLS[index]->Sumw2();
557     hisname.Form("hSumd02SigPt%dLS",i);
558     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
559     fSumd02HistLS[index]->Sumw2();
560     hisname.Form("hSigVertSigPt%dLS",i);
561     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
562     fSigVertHistLS[index]->Sumw2();
563     hisname.Form("hPtMaxSigPt%dLS",i);
564     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
565     fPtMaxHistLS[index]->Sumw2();
566
567     hisname.Form("hDCASigPt%dLS",i);
568     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
569     fDCAHistLS[index]->Sumw2();
570     
571
572
573     index=GetBackgroundHistoIndex(i); 
574     indexLS++;
575     hisname.Form("hBkgPt%d",i);
576     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
577     fMassHist[index]->Sumw2();
578     hisname.Form("hCosPBkgPt%d",i);
579     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
580     fCosPHist[index]->Sumw2();
581     hisname.Form("hDLenBkgPt%d",i);
582     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
583     fDLenHist[index]->Sumw2();
584     hisname.Form("hSumd02BkgPt%d",i);
585     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
586     fSumd02Hist[index]->Sumw2();
587     hisname.Form("hSigVertBkgPt%d",i);
588     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
589     fSigVertHist[index]->Sumw2();
590     hisname.Form("hPtMaxBkgPt%d",i);
591     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
592     fPtMaxHist[index]->Sumw2();
593
594     hisname.Form("hDCABkgPt%d",i);
595     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
596     fDCAHist[index]->Sumw2();
597
598
599     hisname.Form("hBkgPt%dTC",i);
600     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
601     fMassHistTC[index]->Sumw2();
602     hisname.Form("hBkgPt%dTCPlus",i);
603     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
604     fMassHistTCPlus[index]->Sumw2();
605     hisname.Form("hBkgPt%dTCMinus",i);
606     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
607     fMassHistTCMinus[index]->Sumw2();
608
609     hisname.Form("hLSPt%dLCntrip",i);
610     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
611     fMassHistLS[indexLS]->Sumw2();
612     hisname.Form("hLSPt%dTCntrip",i);
613     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
614     fMassHistLSTC[indexLS]->Sumw2();
615
616     
617     hisname.Form("hCosPBkgPt%dLS",i);
618     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
619     fCosPHistLS[index]->Sumw2();
620     hisname.Form("hDLenBkgPt%dLS",i);
621     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
622     fDLenHistLS[index]->Sumw2();
623     hisname.Form("hSumd02BkgPt%dLS",i);
624     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
625     fSumd02HistLS[index]->Sumw2();
626     hisname.Form("hSigVertBkgPt%dLS",i);
627     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
628     fSigVertHistLS[index]->Sumw2();
629     hisname.Form("hPtMaxBkgPt%dLS",i);
630     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
631     fPtMaxHistLS[index]->Sumw2();
632     hisname.Form("hDCABkgPt%dLS",i);
633     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
634     fDCAHistLS[index]->Sumw2();
635     
636
637     indexLS++;
638     hisname.Form("hLSPt%dLCntripsinglecut",i);
639     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
640     fMassHistLS[indexLS]->Sumw2();
641     hisname.Form("hLSPt%dTCntripsinglecut",i);
642     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
643     fMassHistLSTC[indexLS]->Sumw2();
644
645     indexLS++;
646     hisname.Form("hLSPt%dLCspc",i);
647     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
648     fMassHistLS[indexLS]->Sumw2();
649     hisname.Form("hLSPt%dTCspc",i);
650     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
651     fMassHistLSTC[indexLS]->Sumw2();
652   }
653   
654   for(Int_t i=0; i<3*fNPtBins; i++){
655     fOutput->Add(fMassHist[i]);
656      fOutput->Add(fCosPHist[i]);
657     fOutput->Add(fDLenHist[i]);
658     fOutput->Add(fSumd02Hist[i]);
659     fOutput->Add(fSigVertHist[i]);
660     fOutput->Add(fPtMaxHist[i]);
661     fOutput->Add(fDCAHist[i]);
662     fOutput->Add(fMassHistTC[i]);
663     fOutput->Add(fMassHistTCPlus[i]);
664     fOutput->Add(fMassHistTCMinus[i]);
665   }
666   for(Int_t i=0; i<3*fNPtBins&&fDoLS; i++){
667     fOutput->Add(fCosPHistLS[i]);
668     fOutput->Add(fDLenHistLS[i]);
669     fOutput->Add(fSumd02HistLS[i]);
670     fOutput->Add(fSigVertHistLS[i]);
671     fOutput->Add(fPtMaxHistLS[i]);  
672     fOutput->Add(fDCAHistLS[i]);  
673   }
674   for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){
675     fOutput->Add(fMassHistLS[i]);
676     fOutput->Add(fMassHistLSTC[i]);
677   }
678   
679   
680   fHistNEvents = new TH1F("fHistNEvents", "number of events ",9,-0.5,8.5);
681   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
682   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
683   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger");
684   fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events");
685   fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate");
686   fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts");
687   fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts");
688   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
689   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
690
691   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
692   
693   fHistNEvents->Sumw2();
694   fHistNEvents->SetMinimum(0);
695   fOutput->Add(fHistNEvents);
696
697   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
698   fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);  
699   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
700   fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
701   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
702   fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
703
704   fOutput->Add(fPtVsMass);
705   fOutput->Add(fPtVsMassTC);
706   fOutput->Add(fYVsPt);
707   fOutput->Add(fYVsPtTC);
708   fOutput->Add(fYVsPtSig);
709   fOutput->Add(fYVsPtSigTC);
710
711
712   //Counter for Normalization
713   TString normName="NormalizationCounter";
714   AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
715   if(cont)normName=(TString)cont->GetName();
716   fCounter = new AliNormalizationCounter(normName.Data());
717   fCounter->SetRejectPileUp(fRDCutsProduction->GetOptPileUp());
718
719   if(fFillNtuple){
720     OpenFile(4); // 4 is the slot number of the ntuple
721    
722     fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");  
723     
724   }
725   
726   return;
727 }
728
729 //________________________________________________________________________
730 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
731 {
732   // Execute analysis for current event:
733   // heavy flavor candidates association to MC truth
734
735   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
736   
737  
738
739   TClonesArray *array3Prong = 0;
740   TClonesArray *arrayLikeSign =0;
741   if(!aod && AODEvent() && IsStandardAOD()) {
742     // In case there is an AOD handler writing a standard AOD, use the AOD 
743     // event in memory rather than the input (ESD) event.    
744     aod = dynamic_cast<AliAODEvent*> (AODEvent());
745      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
746      // have to taken from the AOD event hold by the AliAODExtension
747     AliAODHandler* aodHandler = (AliAODHandler*) 
748       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
749     if(aodHandler->GetExtensions()) {
750       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
751       AliAODEvent *aodFromExt = ext->GetAOD();
752       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
753       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
754     }
755   } else if(aod) {
756     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
757     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
758   }
759
760   if(!aod || !array3Prong) {
761     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
762     return;
763   }
764   if(!arrayLikeSign && fDoLS) {
765     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
766     return;
767   }
768
769
770   // fix for temporary bug in ESDfilter 
771   // the AODs with null vertex pointer didn't pass the PhysSel
772   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
773   fCounter->StoreEvent(aod,fReadMC);
774   fHistNEvents->Fill(0); // count event
775
776   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
777   Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
778   fHistCentrality[0]->Fill(centrality);
779   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
780   TString trigclass=aod->GetFiredTriggerClasses();
781   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
782   if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3); 
783   if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
784
785   // Post the data already here  
786   PostData(1,fOutput);
787   if(!isEvSel)return;
788
789   fHistCentrality[1]->Fill(centrality);
790   fHistNEvents->Fill(1);
791
792   TClonesArray *arrayMC=0;
793   AliAODMCHeader *mcHeader=0;
794
795   // AOD primary vertex
796   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
797   //    vtx1->Print();
798    TString primTitle = vtx1->GetTitle();
799    //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
800  
801   // load MC particles
802   if(fReadMC){
803     
804     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
805     if(!arrayMC) {
806       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
807       return;
808     }
809     
810   // load MC header
811     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
812     if(!mcHeader) {
813     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
814     return;
815     }
816   }
817   
818   Int_t n3Prong = array3Prong->GetEntriesFast();
819   printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
820   
821   
822   Int_t nOS=0;
823   Int_t index;
824   Int_t pdgDgDplustoKpipi[3]={321,211,211};
825   // 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
826   //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
827   Int_t nSelectedloose=0,nSelectedtight=0;
828   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
829     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
830     fHistNEvents->Fill(4);
831     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
832       fHistNEvents->Fill(8);
833       continue;
834     }
835     Bool_t unsetvtx=kFALSE;
836     if(!d->GetOwnPrimaryVtx()){
837       d->SetOwnPrimaryVtx(vtx1);
838       unsetvtx=kTRUE;
839     }
840
841     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
842
843
844       
845
846       Int_t iPtBin = -1;
847       Double_t ptCand = d->Pt();
848
849       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
850         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
851       }
852       
853       Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
854      
855       
856       Int_t labDp=-1;
857       Float_t deltaPx=0.;
858       Float_t deltaPy=0.;
859       Float_t deltaPz=0.;
860       Float_t truePt=0.;
861       Float_t xDecay=0.;
862       Float_t yDecay=0.;
863       Float_t zDecay=0.;
864       Float_t pdgCode=-2;
865       if(fReadMC){
866         labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
867         if(labDp>=0){
868           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
869           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
870           deltaPx=partDp->Px()-d->Px();
871           deltaPy=partDp->Py()-d->Py();
872           deltaPz=partDp->Pz()-d->Pz();
873           truePt=partDp->Pt();
874           xDecay=dg0->Xv();       
875           yDecay=dg0->Yv();       
876           zDecay=dg0->Zv();
877           pdgCode=TMath::Abs(partDp->GetPdgCode());
878         }else{
879           pdgCode=-1;
880         }
881       }
882       Double_t invMass=d->InvMassDplus();
883       Double_t rapid=d->YDplus();
884       fYVsPt->Fill(ptCand,rapid);
885       if(passTightCuts) fYVsPtTC->Fill(ptCand,rapid);
886       Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
887       if(isFidAcc){
888         fPtVsMass->Fill(invMass,ptCand);
889         if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
890       }
891       Float_t tmp[24];
892       if(fFillNtuple){            
893         tmp[0]=pdgCode;
894         tmp[1]=deltaPx;
895         tmp[2]=deltaPy;
896         tmp[3]=deltaPz;
897         tmp[4]=truePt;
898         tmp[5]=xDecay;    
899         tmp[6]=yDecay;    
900         tmp[7]=zDecay;    
901         tmp[8]=d->PtProng(0);
902         tmp[9]=d->PtProng(1);
903         tmp[10]=d->PtProng(2);
904         tmp[11]=d->Pt();
905         tmp[12]=d->CosPointingAngle();
906         tmp[13]=d->DecayLength();
907         tmp[14]=d->Xv();
908         tmp[15]=d->Yv();
909         tmp[16]=d->Zv();
910         tmp[17]=d->InvMassDplus();
911         tmp[18]=d->GetSigmaVert();
912         tmp[19]=d->Getd0Prong(0);
913         tmp[20]=d->Getd0Prong(1);
914         tmp[21]=d->Getd0Prong(2);
915         tmp[22]=d->GetDCA();
916         tmp[23]=d->Prodd0d0(); 
917         fNtupleDplus->Fill(tmp);
918         PostData(4,fNtupleDplus);
919       }
920       Double_t dlen=d->DecayLength();
921       Double_t cosp=d->CosPointingAngle();
922       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
923       Double_t dca=d->GetDCA();
924       Double_t sigvert=d->GetSigmaVert();         
925       Double_t ptmax=0;
926       for(Int_t i=0;i<3;i++){
927         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
928       }
929       if(iPtBin>=0){
930       
931         index=GetHistoIndex(iPtBin);
932         if(isFidAcc){
933           fHistNEvents->Fill(5);
934           nSelectedloose++;
935           fMassHist[index]->Fill(invMass);
936           fCosPHist[index]->Fill(cosp);
937           fDLenHist[index]->Fill(dlen);
938           fSumd02Hist[index]->Fill(sumD02);
939           fSigVertHist[index]->Fill(sigvert);
940           fPtMaxHist[index]->Fill(ptmax);
941           fDCAHist[index]->Fill(dca);
942           
943           if(passTightCuts){ fHistNEvents->Fill(6);
944             nSelectedtight++;
945             fMassHistTC[index]->Fill(invMass);
946             if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
947             else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
948           }
949         }
950         
951         if(fReadMC){
952           if(labDp>=0) {
953             index=GetSignalHistoIndex(iPtBin);
954             if(isFidAcc){
955               Float_t factor[3]={1.,1.,1.};
956               if(fUseStrangeness){
957                 for(Int_t iprong=0;iprong<3;iprong++){
958                   AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
959                   Int_t labd= trad->GetLabel();
960                   if(labd>=0){
961                     AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
962                     if(dau){
963                       Int_t labm = dau->GetMother();
964                       if(labm>=0){
965                         AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
966                         if(mot){
967                           if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
968                             if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
969                             else factor[iprong]=1./.6;
970                             //          fNentries->Fill(12);
971                           }
972                           if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
973                             factor[iprong]=1./0.25;
974                             //            fNentries->Fill(13);
975                           }//if 3122
976                         }//if(mot)
977                       }//if labm>0
978                     }//if(dau)
979                   }//if labd>=0
980                 }//prong loop
981               }
982               Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
983               fMassHist[index]->Fill(invMass);
984               fCosPHist[index]->Fill(cosp,fact);
985               fDLenHist[index]->Fill(dlen,fact);
986               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];
987               fSumd02Hist[index]->Fill(sumd02s);
988               fSigVertHist[index]->Fill(sigvert,fact);
989               fPtMaxHist[index]->Fill(ptmax,fact);
990               fDCAHist[index]->Fill(dca,fact);
991               if(passTightCuts){
992                 fMassHistTC[index]->Fill(invMass);            
993                 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
994                 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
995               }
996             }       
997             fYVsPtSig->Fill(ptCand,rapid);
998             if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
999           }else{
1000             index=GetBackgroundHistoIndex(iPtBin);
1001             if(isFidAcc){
1002               Float_t factor[3]={1.,1.,1.};
1003               if(fUseStrangeness){
1004                 for(Int_t iprong=0;iprong<3;iprong++){
1005                   AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1006                   Int_t labd= trad->GetLabel();
1007                   if(labd>=0){
1008                     AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1009                     if(dau){
1010                       Int_t labm = dau->GetMother();
1011                       if(labm>=0){
1012                         AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1013                         if(mot){
1014                           if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1015                             if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1016                             else factor[iprong]=1./.6;
1017                             //          fNentries->Fill(12);
1018                           }
1019                           if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1020                             factor[iprong]=1./0.25;
1021                             //            fNentries->Fill(13);
1022                           }//if 3122
1023                         }//if(mot)
1024                       }//if labm>0
1025                     }//if(dau)
1026                   }//if labd>=0
1027                 }//prong loop
1028               }
1029               Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1030               fMassHist[index]->Fill(invMass);
1031               fCosPHist[index]->Fill(cosp,fact);
1032               fDLenHist[index]->Fill(dlen,fact);
1033               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];
1034               fSumd02Hist[index]->Fill(sumd02s);
1035               fSigVertHist[index]->Fill(sigvert,fact);
1036               fPtMaxHist[index]->Fill(ptmax,fact);
1037               fDCAHist[index]->Fill(dca,fact);
1038               if(passTightCuts){
1039                 fMassHistTC[index]->Fill(invMass);
1040                 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1041                 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1042               }
1043             }   
1044           }
1045       
1046           }
1047     
1048     }
1049       
1050     }
1051     if(unsetvtx) d->UnsetOwnPrimaryVtx();
1052   }
1053   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1054   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1055   
1056   //start LS analysis
1057   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1058   
1059   PostData(1,fOutput);    
1060   PostData(3,fCounter);    
1061   return;
1062 }
1063
1064
1065
1066 //________________________________________________________________________
1067 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1068 {
1069   // Terminate analysis
1070   //
1071   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1072
1073   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1074   if (!fOutput) {     
1075     printf("ERROR: fOutput not available\n");
1076     return;
1077   }
1078   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1079
1080   TString hisname;
1081   Int_t index=0;
1082
1083   for(Int_t i=0;i<fNPtBins;i++){
1084     index=GetHistoIndex(i);
1085     
1086     hisname.Form("hMassPt%dTC",i);
1087     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1088   } 
1089     
1090   TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1091   c1->cd();
1092   TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1093   hMassPt->SetLineColor(kBlue);
1094   hMassPt->SetXTitle("M[GeV/c^{2}]"); 
1095   hMassPt->Draw();
1096  
1097   return;
1098 }