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