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