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