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