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