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