a598db2ebe904e9b6370e06004153025087e31b3
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSEDs.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, 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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////
19 //                                                               //
20 //  Analysis task to produce Ds candidates mass spectra          //
21 // Origin: F.Prino, Torino, prino@to.infn.it                     //
22 //                                                               //
23 ///////////////////////////////////////////////////////////////////
24
25 #include <TClonesArray.h>
26 #include <TNtuple.h>
27 #include <TList.h>
28 #include <TString.h>
29 #include <TH1F.h>
30 #include <TMath.h>
31 #include <TDatabasePDG.h>
32
33 #include "AliAnalysisManager.h"
34 #include "AliAODHandler.h"
35 #include "AliAODEvent.h"
36 #include "AliAODVertex.h"
37 #include "AliAODTrack.h"
38 #include "AliAODMCHeader.h"
39 #include "AliAODMCParticle.h"
40 #include "AliAODRecoDecayHF3Prong.h"
41 #include "AliAnalysisVertexingHF.h"
42 #include "AliRDHFCutsDstoKKpi.h"
43 #include "AliAnalysisTaskSE.h"
44 #include "AliAnalysisTaskSEDs.h"
45 #include "AliNormalizationCounter.h"
46
47 ClassImp(AliAnalysisTaskSEDs)
48
49
50 //________________________________________________________________________
51 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
52   AliAnalysisTaskSE(),
53   fOutput(0), 
54   fHistNEvents(0),
55   fPtVsMass(0),
56   fYVsPt(0),
57   fYVsPtSig(0),
58   fNtupleDs(0),
59   fFillNtuple(0),
60   fReadMC(kFALSE),
61   fDoCutVarHistos(kTRUE),
62   fUseSelectionBit(kFALSE),
63   fNPtBins(0),
64   fListCuts(0),
65   fMassRange(0.4),
66   fMassBinSize(0.002),
67   fCounter(0),
68   fProdCuts(0),
69   fAnalysisCuts(0)
70 {
71   // Default constructor
72 }
73
74 //________________________________________________________________________
75 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name, AliRDHFCutsDstoKKpi* productioncuts, AliRDHFCutsDstoKKpi* analysiscuts,Int_t fillNtuple):
76   AliAnalysisTaskSE(name),
77   fOutput(0),
78   fHistNEvents(0),
79   fPtVsMass(0),
80   fYVsPt(0),
81   fYVsPtSig(0),
82   fNtupleDs(0),
83   fFillNtuple(fillNtuple),
84   fReadMC(kFALSE),
85   fDoCutVarHistos(kTRUE),
86   fUseSelectionBit(kFALSE),
87   fNPtBins(0),
88   fListCuts(0),
89   fMassRange(0.4),
90   fMassBinSize(0.002),
91   fCounter(0),
92   fProdCuts(productioncuts),
93   fAnalysisCuts(analysiscuts)
94 {
95   // Default constructor
96   // Output slot #1 writes into a TList container
97   Int_t nptbins=fAnalysisCuts->GetNPtBins();
98   Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
99   SetPtBins(nptbins,ptlim);
100  
101   DefineOutput(1,TList::Class());  //My private output
102
103   DefineOutput(2,TList::Class());
104
105   DefineOutput(3,AliNormalizationCounter::Class());
106   
107   if(fFillNtuple>0){
108     // Output slot #4 writes into a TNtuple container
109     DefineOutput(4,TNtuple::Class());  //My private output
110   }
111   
112 }
113
114 //________________________________________________________________________
115 void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
116   // define pt bins for analysis
117   if(n>kMaxPtBins){
118     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
119     fNPtBins=kMaxPtBins;
120     fPtLimits[0]=0.;
121     fPtLimits[1]=1.;
122     fPtLimits[2]=3.;
123     fPtLimits[3]=5.;
124     fPtLimits[4]=10.;
125     for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
126   }else{
127     fNPtBins=n;
128     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
129     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
130   }
131   if(fDebug > 1){
132     printf("Number of Pt bins = %d\n",fNPtBins);
133     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);    
134   }
135 }
136 //________________________________________________________________________
137 AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
138 {
139   // Destructor
140   if (fOutput) {
141     delete fOutput;
142     fOutput = 0;
143   }
144
145   if(fHistNEvents){
146     delete fHistNEvents;
147     fHistNEvents=0;
148   } 
149  
150   if (fListCuts) {
151     delete fListCuts;
152     fListCuts = 0;
153   }
154
155   for(Int_t i=0;i<4*fNPtBins;i++){
156     
157     if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
158     if(fMassHistPhi[i]){ delete fMassHistPhi[i]; fMassHistPhi[i]=0;}
159     if(fMassHistK0st[i]){ delete fMassHistK0st[i]; fMassHistK0st[i]=0;}
160     if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
161     if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
162     if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
163     if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
164     if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
165     if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
166     if(fPtProng0Hist[i]){ delete fPtProng0Hist[i]; fPtProng0Hist[i]=0;}
167     if(fPtProng1Hist[i]){ delete fPtProng1Hist[i]; fPtProng1Hist[i]=0;}
168     if(fPtProng2Hist[i]){ delete fPtProng2Hist[i]; fPtProng2Hist[i]=0;}
169     if(fDalitz[i]){ delete fDalitz[i]; fDalitz[i]=0;}
170     if(fDalitzPhi[i]){ delete fDalitzPhi[i]; fDalitzPhi[i]=0;}
171     if(fDalitzK0st[i]){ delete fDalitzK0st[i]; fDalitzK0st[i]=0;}
172
173   }
174
175   if(fPtVsMass){
176     delete fPtVsMass;
177     fPtVsMass=0;
178   }
179   if(fYVsPt){
180     delete fYVsPt;
181     fYVsPt=0;
182   }
183   if(fYVsPtSig){
184     delete fYVsPtSig;
185     fYVsPtSig=0;
186   }
187   if(fNtupleDs){
188     delete fNtupleDs;
189     fNtupleDs=0;
190   }
191   if(fCounter){
192     delete fCounter;
193     fCounter = 0;
194   }
195   
196   if (fProdCuts) {
197     delete fProdCuts;
198     fProdCuts = 0;
199   }
200   if (fAnalysisCuts) {
201     delete fAnalysisCuts;
202     fAnalysisCuts = 0;
203   }
204 }  
205
206 //________________________________________________________________________
207 void AliAnalysisTaskSEDs::Init()
208 {
209   // Initialization
210
211   if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
212
213   fListCuts=new TList();
214   fListCuts->SetOwner();
215   fListCuts->SetName("CutObjects");
216
217   AliRDHFCutsDstoKKpi *production = new AliRDHFCutsDstoKKpi(*fProdCuts);
218   production->SetName("ProductionCuts");
219   AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi(*fAnalysisCuts);
220   analysis->SetName("AnalysisCuts");
221   
222   fListCuts->Add(production);
223   fListCuts->Add(analysis);
224   PostData(2,fListCuts);
225   return;
226 }
227
228 //________________________________________________________________________
229 void AliAnalysisTaskSEDs::UserCreateOutputObjects()
230 {
231   // Create the output container
232   //
233   if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
234
235   // Several histograms are more conveniently managed in a TList
236   fOutput = new TList();
237   fOutput->SetOwner();
238   fOutput->SetName("OutputHistos");
239
240   Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
241   
242   Int_t nInvMassBins=(Int_t)(fMassRange/fMassBinSize+0.5);
243   if(nInvMassBins%2==1) nInvMassBins++;
244   // Double_t minMass=massDs-1.0*nInvMassBins*fMassBinSize;
245   Double_t minMass=massDs-0.5*nInvMassBins*fMassBinSize;
246   //  Double_t maxMass=massDs+1.0*nInvMassBins*fMassBinSize;
247   Double_t maxMass=massDs+0.5*nInvMassBins*fMassBinSize;
248
249   TString hisname;
250   TString htype;
251   Int_t index;
252   for(Int_t iType=0; iType<4; iType++){
253     for(Int_t i=0;i<fNPtBins;i++){
254       if(iType==0){
255         htype="All";
256         index=GetHistoIndex(i);
257       }else if(iType==1){ 
258         htype="Sig";
259         index=GetSignalHistoIndex(i);
260       }else if(iType==2){ 
261         htype="Bkg";
262         index=GetBackgroundHistoIndex(i);
263       }else{ 
264         htype="ReflSig";
265         index=GetReflSignalHistoIndex(i);
266       }
267       hisname.Form("hMass%sPt%d",htype.Data(),i);
268       fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
269       fMassHist[index]->Sumw2();
270       hisname.Form("hMass%sPt%dphi",htype.Data(),i);
271       fMassHistPhi[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
272       fMassHistPhi[index]->Sumw2();
273       hisname.Form("hMass%sPt%dk0st",htype.Data(),i);
274       fMassHistK0st[index]=new TH1F(hisname.Data(),hisname.Data(),nInvMassBins,minMass,maxMass);
275       fMassHistK0st[index]->Sumw2();
276       hisname.Form("hCosP%sPt%d",htype.Data(),i);
277       fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
278       fCosPHist[index]->Sumw2();
279       hisname.Form("hDLen%sPt%d",htype.Data(),i);
280       fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
281       fDLenHist[index]->Sumw2();
282       hisname.Form("hSumd02%sPt%d",htype.Data(),i);
283       fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
284       fSumd02Hist[index]->Sumw2();
285       hisname.Form("hSigVert%sPt%d",htype.Data(),i);
286       fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
287       fSigVertHist[index]->Sumw2();
288       hisname.Form("hPtMax%sPt%d",htype.Data(),i);
289       fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
290       fPtMaxHist[index]->Sumw2();
291       hisname.Form("hPtCand%sPt%d",htype.Data(),i);
292       fPtCandHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,20.);
293       fPtCandHist[index]->Sumw2();
294       hisname.Form("hDCA%sPt%d",htype.Data(),i);
295       fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
296       fDCAHist[index]->Sumw2();
297       hisname.Form("hPtProng0%sPt%d",htype.Data(),i);
298       fPtProng0Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
299       fPtProng0Hist[index]->Sumw2();
300       hisname.Form("hPtProng1%sPt%d",htype.Data(),i);
301       fPtProng1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
302       fPtProng1Hist[index]->Sumw2();
303       hisname.Form("hPtProng2%sPt%d",htype.Data(),i);
304       fPtProng2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.0,20.);
305       fPtProng2Hist[index]->Sumw2();
306       hisname.Form("hDalitz%sPt%d",htype.Data(),i);
307       fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
308       fDalitz[index]->Sumw2();
309       hisname.Form("hDalitz%sPt%dphi",htype.Data(),i);
310       fDalitzPhi[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
311       fDalitzPhi[index]->Sumw2();
312       hisname.Form("hDalitz%sPt%dk0st",htype.Data(),i);
313       fDalitzK0st[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
314       fDalitzK0st[index]->Sumw2();
315     }
316   }
317
318   for(Int_t i=0; i<4*fNPtBins; i++){
319     fOutput->Add(fMassHist[i]);
320     fOutput->Add(fMassHistPhi[i]);
321     fOutput->Add(fMassHistK0st[i]);
322     fOutput->Add(fCosPHist[i]);
323     fOutput->Add(fDLenHist[i]);
324     fOutput->Add(fSumd02Hist[i]);
325     fOutput->Add(fSigVertHist[i]);
326     fOutput->Add(fPtMaxHist[i]);
327     fOutput->Add(fPtCandHist[i]);
328     fOutput->Add(fDCAHist[i]);
329     fOutput->Add(fPtProng0Hist[i]);
330     fOutput->Add(fPtProng1Hist[i]);
331     fOutput->Add(fPtProng2Hist[i]);
332     fOutput->Add(fDalitz[i]);
333     fOutput->Add(fDalitzPhi[i]);
334     fOutput->Add(fDalitzK0st[i]);
335   }
336
337   fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",64,-0.5,63.5);
338   fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",64,-0.5,63.5);
339   fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",64,-0.5,63.5);
340   fChanHist[3] = new TH1F("hChanReflSig", "KKpi and piKK candidates",64,-0.5,63.5);
341   for(Int_t i=0;i<4;i++){
342     fChanHist[i]->Sumw2();
343     fChanHist[i]->SetMinimum(0);
344     fOutput->Add(fChanHist[i]);
345   }
346
347   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
348   fHistNEvents->Sumw2();
349   fHistNEvents->SetMinimum(0);
350   fOutput->Add(fHistNEvents);
351
352   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
353   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
354   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
355
356   fOutput->Add(fPtVsMass);
357   fOutput->Add(fYVsPt);
358   fOutput->Add(fYVsPtSig);
359
360   //Counter for Normalization
361   fCounter = new AliNormalizationCounter("NormalizationCounter");
362
363   PostData(1,fOutput); 
364   PostData(3,fCounter);   
365   
366   if(fFillNtuple>0){
367     OpenFile(4); // 4 is the slot number of the ntuple
368     
369     fNtupleDs = new TNtuple("fNtupleDs","Ds","labDs:retcode:pdgcode0:Pt0:Pt1:Pt2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMassKKpi:InvMasspiKK:sigvert:d00:d01:d02:dca:d0square:InvMassPhiKKpi:InvMassPhipiKK:cosinePiDsFrameKKpi:cosinePiDsFramepiKK:cosineKPhiFrameKKpi:cosineKPhiFramepiKK"); 
370     
371   }
372   
373
374   return;
375 }
376
377 //________________________________________________________________________
378 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
379 {
380   // Ds selection for current event, fill mass histos and selecetion variable histo
381   // separate signal and backgound if fReadMC is activated
382
383   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
384
385   TClonesArray *array3Prong = 0;
386   if(!aod && AODEvent() && IsStandardAOD()) {
387     // In case there is an AOD handler writing a standard AOD, use the AOD 
388     // event in memory rather than the input (ESD) event.    
389     aod = dynamic_cast<AliAODEvent*> (AODEvent());
390     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
391     // have to taken from the AOD event hold by the AliAODExtension
392     AliAODHandler* aodHandler = (AliAODHandler*) 
393       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
394     if(aodHandler->GetExtensions()) {
395       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
396       AliAODEvent *aodFromExt = ext->GetAOD();
397       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
398     }
399   } else if(aod) {
400     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
401   }
402
403   if(!aod || !array3Prong) {
404     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
405     return;
406   }
407
408
409   // fix for temporary bug in ESDfilter 
410   // the AODs with null vertex pointer didn't pass the PhysSel
411   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
412
413
414   fHistNEvents->Fill(0); // count event
415   // Post the data already here
416   PostData(1,fOutput);
417
418   fCounter->StoreEvent(aod,fProdCuts,fReadMC);
419   //fCounter->StoreEvent(aod,fReadMC);
420
421   TClonesArray *arrayMC=0;
422   AliAODMCHeader *mcHeader=0;
423
424   // AOD primary vertex
425   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
426   //    vtx1->Print();
427   
428   // load MC particles
429   if(fReadMC){
430     
431     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
432     if(!arrayMC) {
433       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
434       return;
435     }
436     
437   // load MC header
438     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
439     if(!mcHeader) {
440       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
441       return;
442     }
443   }
444   
445   Int_t n3Prong = array3Prong->GetEntriesFast();
446   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
447   
448   
449   Int_t pdgDstoKKpi[3]={321,321,211};
450   Int_t nSelectedloose=0;
451   Int_t nSelectedtight=0;
452
453   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
454   
455     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
456     
457     if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))) continue;
458     
459     Bool_t unsetvtx=kFALSE;
460     if(!d->GetOwnPrimaryVtx()){
461       d->SetOwnPrimaryVtx(vtx1);
462       unsetvtx=kTRUE;
463     }
464
465     Double_t ptCand = d->Pt();
466     Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
467     Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
468     Double_t rapid=d->YDs(); 
469     fYVsPt->Fill(ptCand,rapid);
470
471     Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
472     if(isFidAcc){
473       nSelectedloose++;
474       if(retCodeAnalysisCuts>0)nSelectedtight++;
475     }
476     if(retCodeAnalysisCuts<=0) continue;
477     if(!isFidAcc) continue;
478
479     Int_t index=GetHistoIndex(iPtBin);
480     fPtCandHist[index]->Fill(ptCand);
481
482     Int_t isKKpi=retCodeAnalysisCuts&1;
483     Int_t ispiKK=retCodeAnalysisCuts&2;
484     Int_t isPhiKKpi=retCodeAnalysisCuts&4;
485     Int_t isPhipiKK=retCodeAnalysisCuts&8;
486     Int_t isK0starKKpi=retCodeAnalysisCuts&16;    
487     Int_t isK0starpiKK=retCodeAnalysisCuts&32;
488
489     fChanHist[0]->Fill(retCodeAnalysisCuts);
490  
491     Int_t indexMCKKpi=-1;
492     Int_t indexMCpiKK=-1;
493     Int_t labDs=-1;
494     Int_t pdgCode0=-999;
495     
496     if(fReadMC){
497       labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
498       if(labDs>=0){
499         Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
500         AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
501         pdgCode0=TMath::Abs(p->GetPdgCode());
502
503         if(isKKpi){
504           if(pdgCode0==321) {     
505             indexMCKKpi=GetSignalHistoIndex(iPtBin);
506             fYVsPtSig->Fill(ptCand,rapid);
507             fChanHist[1]->Fill(retCodeAnalysisCuts);
508           }else{
509             indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
510             fChanHist[3]->Fill(retCodeAnalysisCuts);
511           }
512         }
513         if(ispiKK){
514           if(pdgCode0==211) {     
515             indexMCpiKK=GetSignalHistoIndex(iPtBin);
516             fYVsPtSig->Fill(ptCand,rapid);
517             fChanHist[1]->Fill(retCodeAnalysisCuts);
518           }else{
519             indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
520             fChanHist[3]->Fill(retCodeAnalysisCuts);
521           }
522         }
523       }else{
524         indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
525         indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
526         fChanHist[2]->Fill(retCodeAnalysisCuts);
527       }
528     }
529
530     if(isKKpi){
531       Double_t invMass=d->InvMassDsKKpi();
532       fMassHist[index]->Fill(invMass);
533       fPtVsMass->Fill(invMass,ptCand);
534       if(isPhiKKpi) fMassHistPhi[index]->Fill(invMass); 
535       if(isK0starKKpi) fMassHistK0st[index]->Fill(invMass);
536       if(fReadMC  && indexMCKKpi!=-1){
537         fMassHist[indexMCKKpi]->Fill(invMass);
538         if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass);
539         if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass);       
540       }
541     }
542     if(ispiKK){
543       Double_t invMass=d->InvMassDspiKK();
544       fMassHist[index]->Fill(invMass);
545       fPtVsMass->Fill(invMass,ptCand);
546       if(isPhipiKK) fMassHistPhi[index]->Fill(invMass);
547       if(isK0starpiKK) fMassHistK0st[index]->Fill(invMass);
548       if(fReadMC  && indexMCpiKK!=-1){
549         fMassHist[indexMCpiKK]->Fill(invMass);
550         if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass);
551         if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass);      
552       }
553     }
554
555     if(fDoCutVarHistos){
556       Double_t dlen=d->DecayLength();
557       Double_t cosp=d->CosPointingAngle();
558       Double_t pt0=d->PtProng(0);
559       Double_t pt1=d->PtProng(1);
560       Double_t pt2=d->PtProng(2);
561       Double_t sigvert=d->GetSigmaVert();
562       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
563       Double_t dca=d->GetDCA();      
564       Double_t ptmax=0;
565       for(Int_t i=0;i<3;i++){
566         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
567       }
568       fCosPHist[index]->Fill(cosp);
569       fDLenHist[index]->Fill(dlen);
570       fSigVertHist[index]->Fill(sigvert);
571       fSumd02Hist[index]->Fill(sumD02);
572       fPtMaxHist[index]->Fill(ptmax);
573       fDCAHist[index]->Fill(dca);
574       fPtProng0Hist[index]->Fill(pt0);
575       fPtProng1Hist[index]->Fill(pt1);
576       fPtProng2Hist[index]->Fill(pt2);
577       if(isKKpi){
578         Double_t massKK=d->InvMass2Prongs(0,1,321,321);
579         Double_t massKp=d->InvMass2Prongs(1,2,321,211);
580         fDalitz[index]->Fill(massKK,massKp);
581         if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
582         if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
583         if(fReadMC && indexMCKKpi!=-1){
584           fDalitz[indexMCKKpi]->Fill(massKK,massKp);
585           if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
586           if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
587           fCosPHist[indexMCKKpi]->Fill(cosp);
588           fDLenHist[indexMCKKpi]->Fill(dlen);
589           fSigVertHist[indexMCKKpi]->Fill(sigvert);
590           fSumd02Hist[indexMCKKpi]->Fill(sumD02);
591           fPtMaxHist[indexMCKKpi]->Fill(ptmax);
592           fPtCandHist[indexMCKKpi]->Fill(ptCand);
593           fDCAHist[indexMCKKpi]->Fill(dca);
594           fPtProng0Hist[indexMCKKpi]->Fill(pt0);
595           fPtProng1Hist[indexMCKKpi]->Fill(pt1);
596           fPtProng2Hist[indexMCKKpi]->Fill(pt2);
597         }       
598       }
599       if(ispiKK){
600         Double_t massKK=d->InvMass2Prongs(1,2,321,321);
601         Double_t massKp=d->InvMass2Prongs(0,1,211,321);
602         fDalitz[index]->Fill(massKK,massKp);
603         if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
604         if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
605
606         
607         if(fReadMC && indexMCpiKK!=-1){
608           fDalitz[indexMCpiKK]->Fill(massKK,massKp);
609           if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
610           if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
611           fCosPHist[indexMCpiKK]->Fill(cosp);
612           fDLenHist[indexMCpiKK]->Fill(dlen);
613           fSigVertHist[indexMCpiKK]->Fill(sigvert);
614           fSumd02Hist[indexMCpiKK]->Fill(sumD02);
615           fPtMaxHist[indexMCpiKK]->Fill(ptmax);
616           fPtCandHist[indexMCpiKK]->Fill(ptCand);
617           fDCAHist[indexMCpiKK]->Fill(dca);
618           fPtProng0Hist[indexMCpiKK]->Fill(pt0);
619           fPtProng1Hist[indexMCpiKK]->Fill(pt1);
620           fPtProng2Hist[indexMCpiKK]->Fill(pt2);
621         }
622       }
623      
624    
625     }
626     
627     Float_t tmp[26];
628     if(fFillNtuple>0){
629       
630       if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
631         
632         tmp[0]=Float_t(labDs);
633         tmp[1]=Float_t(retCodeAnalysisCuts);
634         tmp[2]=Float_t(pdgCode0);  
635         tmp[3]=d->PtProng(0);
636         tmp[4]=d->PtProng(1);
637         tmp[5]=d->PtProng(2);
638         tmp[6]=d->Pt();
639         tmp[7]=d->CosPointingAngle();
640         tmp[8]=d->DecayLength();
641         tmp[9]=d->Xv();
642         tmp[10]=d->Yv();
643         tmp[11]=d->Zv();
644         tmp[12]=d->InvMassDsKKpi();
645         tmp[13]=d->InvMassDspiKK();
646         tmp[14]=d->GetSigmaVert();
647         tmp[15]=d->Getd0Prong(0);
648         tmp[16]=d->Getd0Prong(1);
649         tmp[17]=d->Getd0Prong(2);
650         tmp[18]=d->GetDCA();
651         tmp[19]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
652         tmp[20]=d->InvMass2Prongs(0,1,321,321);
653         tmp[21]=d->InvMass2Prongs(1,2,321,321);
654         tmp[22]=d->CosPiDsLabFrameKKpi();          
655         tmp[23]=d->CosPiDsLabFramepiKK();       
656         tmp[24]=d->CosPiKPhiRFrameKKpi();          
657         tmp[25]=d->CosPiKPhiRFramepiKK();       
658         
659         
660         fNtupleDs->Fill(tmp);
661         PostData(4,fNtupleDs);
662       }  
663     }
664     
665     if(unsetvtx) d->UnsetOwnPrimaryVtx();
666   }
667  
668   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
669   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
670
671   PostData(1,fOutput); 
672   PostData(3,fCounter);    
673
674   return;
675 }
676
677 //_________________________________________________________________
678
679 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
680 {
681   // Terminate analysis
682   //
683   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
684   fOutput = dynamic_cast<TList*> (GetOutputData(1));
685   if (!fOutput) {     
686     printf("ERROR: fOutput not available\n");
687     return;
688   }
689   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
690   if(fHistNEvents){
691     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
692   }else{
693     printf("ERROR: fHistNEvents not available\n");
694     return;
695   }
696   return;
697 }
698