Add recalculation of primary vertex without daughters, updated histogram with event...
[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 events ",12,-0.5,11.5);
348   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
349   fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
350   fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
351   fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
352   fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
353   fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
354   fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
355   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
356   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of candidate");
357   fHistNEvents->GetXaxis()->SetBinLabel(10,"no. of Ds after loose cuts");
358   fHistNEvents->GetXaxis()->SetBinLabel(11,"no. of Ds after tight cuts");
359   fHistNEvents->GetXaxis()->SetBinLabel(12,"no. of cand wo bitmask");
360
361   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
362
363   fHistNEvents->Sumw2();
364   fHistNEvents->SetMinimum(0);
365   fOutput->Add(fHistNEvents);
366
367   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nInvMassBins,minMass,maxMass,40,0.,20.);
368   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
369   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
370
371   fOutput->Add(fPtVsMass);
372   fOutput->Add(fYVsPt);
373   fOutput->Add(fYVsPtSig);
374
375   //Counter for Normalization
376   fCounter = new AliNormalizationCounter("NormalizationCounter");
377
378   PostData(1,fOutput); 
379   PostData(3,fCounter);   
380   
381   if(fFillNtuple>0){
382     OpenFile(4); // 4 is the slot number of the ntuple
383     
384     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"); 
385     
386   }
387   
388
389   return;
390 }
391
392 //________________________________________________________________________
393 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
394 {
395   // Ds selection for current event, fill mass histos and selecetion variable histo
396   // separate signal and backgound if fReadMC is activated
397
398   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
399
400   TClonesArray *array3Prong = 0;
401   if(!aod && AODEvent() && IsStandardAOD()) {
402     // In case there is an AOD handler writing a standard AOD, use the AOD 
403     // event in memory rather than the input (ESD) event.    
404     aod = dynamic_cast<AliAODEvent*> (AODEvent());
405     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
406     // have to taken from the AOD event hold by the AliAODExtension
407     AliAODHandler* aodHandler = (AliAODHandler*) 
408       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
409     if(aodHandler->GetExtensions()) {
410       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
411       AliAODEvent *aodFromExt = ext->GetAOD();
412       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
413     }
414   } else if(aod) {
415     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
416   }
417
418   if(!aod || !array3Prong) {
419     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
420     return;
421   }
422
423
424   // fix for temporary bug in ESDfilter 
425   // the AODs with null vertex pointer didn't pass the PhysSel
426   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
427
428
429   fHistNEvents->Fill(0); // count event
430   // Post the data already here
431   PostData(1,fOutput);
432   
433   fCounter->StoreEvent(aod,fProdCuts,fReadMC);
434   
435
436   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
437   if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
438   if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
439   if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
440   if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
441   if(fAnalysisCuts->IsEventRejectedDueToPileupSPD())fHistNEvents->Fill(6);
442   if(fAnalysisCuts->IsEventRejectedDueToCentrality())fHistNEvents->Fill(7);
443   
444   
445   
446   if(!isEvSel)return;
447   
448   fHistNEvents->Fill(1);
449
450   TClonesArray *arrayMC=0;
451   AliAODMCHeader *mcHeader=0;
452
453   // AOD primary vertex
454   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
455   //    vtx1->Print();
456   
457   // load MC particles
458   if(fReadMC){
459     
460     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
461     if(!arrayMC) {
462       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
463       return;
464     }
465     
466   // load MC header
467     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
468     if(!mcHeader) {
469       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
470       return;
471     }
472   }
473   
474   Int_t n3Prong = array3Prong->GetEntriesFast();
475   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
476   
477   
478   Int_t pdgDstoKKpi[3]={321,321,211};
479   Int_t nSelectedloose=0;
480   Int_t nSelectedtight=0;
481
482   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
483   
484     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
485     fHistNEvents->Fill(8);
486     
487     if(fUseSelectionBit && !(d->HasSelectionBit(AliRDHFCuts::kDsCuts))){
488       fHistNEvents->Fill(11);
489       continue;
490     }
491     
492     Bool_t unsetvtx=kFALSE;
493     if(!d->GetOwnPrimaryVtx()){
494       d->SetOwnPrimaryVtx(vtx1);
495       unsetvtx=kTRUE;
496     }
497     
498     Bool_t recVtx=kFALSE;
499     AliAODVertex *origownvtx=0x0;
500     Int_t retCodeProdCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
501    
502     if(retCodeProdCuts) {
503       if(fProdCuts->GetIsPrimaryWithoutDaughters()){
504             if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());      
505             if(fProdCuts->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
506             else fProdCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
507       }
508     }  
509     
510     Double_t ptCand = d->Pt();
511     Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
512     Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
513     Double_t rapid=d->YDs(); 
514     fYVsPt->Fill(ptCand,rapid);
515
516     Bool_t isFidAcc=fAnalysisCuts->IsInFiducialAcceptance(ptCand,rapid);
517     
518     if(retCodeProdCuts>0){
519       if(isFidAcc){
520         nSelectedloose++;
521         fHistNEvents->Fill(9);
522         if(retCodeAnalysisCuts>0)nSelectedtight++;
523       }
524     }
525   
526     if(retCodeAnalysisCuts<=0) continue;
527     if(!isFidAcc) continue;
528     fHistNEvents->Fill(10);
529     
530     Int_t index=GetHistoIndex(iPtBin);
531     fPtCandHist[index]->Fill(ptCand);
532
533     Int_t isKKpi=retCodeAnalysisCuts&1;
534     Int_t ispiKK=retCodeAnalysisCuts&2;
535     Int_t isPhiKKpi=retCodeAnalysisCuts&4;
536     Int_t isPhipiKK=retCodeAnalysisCuts&8;
537     Int_t isK0starKKpi=retCodeAnalysisCuts&16;    
538     Int_t isK0starpiKK=retCodeAnalysisCuts&32;
539
540     fChanHist[0]->Fill(retCodeAnalysisCuts);
541  
542     Int_t indexMCKKpi=-1;
543     Int_t indexMCpiKK=-1;
544     Int_t labDs=-1;
545     Int_t pdgCode0=-999;
546     
547     if(fReadMC){
548       labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
549       if(labDs>=0){
550         Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
551         AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(TMath::Abs(labDau0));
552         pdgCode0=TMath::Abs(p->GetPdgCode());
553
554         if(isKKpi){
555           if(pdgCode0==321) {     
556             indexMCKKpi=GetSignalHistoIndex(iPtBin);
557             fYVsPtSig->Fill(ptCand,rapid);
558             fChanHist[1]->Fill(retCodeAnalysisCuts);
559           }else{
560             indexMCKKpi=GetReflSignalHistoIndex(iPtBin);
561             fChanHist[3]->Fill(retCodeAnalysisCuts);
562           }
563         }
564         if(ispiKK){
565           if(pdgCode0==211) {     
566             indexMCpiKK=GetSignalHistoIndex(iPtBin);
567             fYVsPtSig->Fill(ptCand,rapid);
568             fChanHist[1]->Fill(retCodeAnalysisCuts);
569           }else{
570             indexMCpiKK=GetReflSignalHistoIndex(iPtBin);
571             fChanHist[3]->Fill(retCodeAnalysisCuts);
572           }
573         }
574       }else{
575         indexMCpiKK=GetBackgroundHistoIndex(iPtBin);
576         indexMCKKpi=GetBackgroundHistoIndex(iPtBin);
577         fChanHist[2]->Fill(retCodeAnalysisCuts);
578       }
579     }
580
581     if(isKKpi){
582       Double_t invMass=d->InvMassDsKKpi();
583       fMassHist[index]->Fill(invMass);
584       fPtVsMass->Fill(invMass,ptCand);
585       if(isPhiKKpi) fMassHistPhi[index]->Fill(invMass); 
586       if(isK0starKKpi) fMassHistK0st[index]->Fill(invMass);
587       if(fReadMC  && indexMCKKpi!=-1){
588         fMassHist[indexMCKKpi]->Fill(invMass);
589         if(isPhiKKpi) fMassHistPhi[indexMCKKpi]->Fill(invMass);
590         if(isK0starKKpi) fMassHistK0st[indexMCKKpi]->Fill(invMass);       
591       }
592     }
593     if(ispiKK){
594       Double_t invMass=d->InvMassDspiKK();
595       fMassHist[index]->Fill(invMass);
596       fPtVsMass->Fill(invMass,ptCand);
597       if(isPhipiKK) fMassHistPhi[index]->Fill(invMass);
598       if(isK0starpiKK) fMassHistK0st[index]->Fill(invMass);
599       if(fReadMC  && indexMCpiKK!=-1){
600         fMassHist[indexMCpiKK]->Fill(invMass);
601         if(isPhipiKK) fMassHistPhi[indexMCpiKK]->Fill(invMass);
602         if(isK0starpiKK) fMassHistK0st[indexMCpiKK]->Fill(invMass);      
603       }
604     }
605
606     if(fDoCutVarHistos){
607       Double_t dlen=d->DecayLength();
608       Double_t cosp=d->CosPointingAngle();
609       Double_t pt0=d->PtProng(0);
610       Double_t pt1=d->PtProng(1);
611       Double_t pt2=d->PtProng(2);
612       Double_t sigvert=d->GetSigmaVert();
613       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
614       Double_t dca=d->GetDCA();      
615       Double_t ptmax=0;
616       for(Int_t i=0;i<3;i++){
617         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
618       }
619       fCosPHist[index]->Fill(cosp);
620       fDLenHist[index]->Fill(dlen);
621       fSigVertHist[index]->Fill(sigvert);
622       fSumd02Hist[index]->Fill(sumD02);
623       fPtMaxHist[index]->Fill(ptmax);
624       fDCAHist[index]->Fill(dca);
625       fPtProng0Hist[index]->Fill(pt0);
626       fPtProng1Hist[index]->Fill(pt1);
627       fPtProng2Hist[index]->Fill(pt2);
628       if(isKKpi){
629         Double_t massKK=d->InvMass2Prongs(0,1,321,321);
630         Double_t massKp=d->InvMass2Prongs(1,2,321,211);
631         fDalitz[index]->Fill(massKK,massKp);
632         if(isPhiKKpi) fDalitzPhi[index]->Fill(massKK,massKp);
633         if(isK0starKKpi) fDalitzK0st[index]->Fill(massKK,massKp);
634         if(fReadMC && indexMCKKpi!=-1){
635           fDalitz[indexMCKKpi]->Fill(massKK,massKp);
636           if(isPhiKKpi) fDalitzPhi[indexMCKKpi]->Fill(massKK,massKp);
637           if(isK0starKKpi) fDalitzK0st[indexMCKKpi]->Fill(massKK,massKp);
638           fCosPHist[indexMCKKpi]->Fill(cosp);
639           fDLenHist[indexMCKKpi]->Fill(dlen);
640           fSigVertHist[indexMCKKpi]->Fill(sigvert);
641           fSumd02Hist[indexMCKKpi]->Fill(sumD02);
642           fPtMaxHist[indexMCKKpi]->Fill(ptmax);
643           fPtCandHist[indexMCKKpi]->Fill(ptCand);
644           fDCAHist[indexMCKKpi]->Fill(dca);
645           fPtProng0Hist[indexMCKKpi]->Fill(pt0);
646           fPtProng1Hist[indexMCKKpi]->Fill(pt1);
647           fPtProng2Hist[indexMCKKpi]->Fill(pt2);
648         }       
649       }
650       if(ispiKK){
651         Double_t massKK=d->InvMass2Prongs(1,2,321,321);
652         Double_t massKp=d->InvMass2Prongs(0,1,211,321);
653         fDalitz[index]->Fill(massKK,massKp);
654         if(isPhipiKK) fDalitzPhi[index]->Fill(massKK,massKp);
655         if(isK0starpiKK) fDalitzK0st[index]->Fill(massKK,massKp);
656
657         
658         if(fReadMC && indexMCpiKK!=-1){
659           fDalitz[indexMCpiKK]->Fill(massKK,massKp);
660           if(isPhipiKK) fDalitzPhi[indexMCpiKK]->Fill(massKK,massKp);
661           if(isK0starpiKK) fDalitzK0st[indexMCpiKK]->Fill(massKK,massKp);
662           fCosPHist[indexMCpiKK]->Fill(cosp);
663           fDLenHist[indexMCpiKK]->Fill(dlen);
664           fSigVertHist[indexMCpiKK]->Fill(sigvert);
665           fSumd02Hist[indexMCpiKK]->Fill(sumD02);
666           fPtMaxHist[indexMCpiKK]->Fill(ptmax);
667           fPtCandHist[indexMCpiKK]->Fill(ptCand);
668           fDCAHist[indexMCpiKK]->Fill(dca);
669           fPtProng0Hist[indexMCpiKK]->Fill(pt0);
670           fPtProng1Hist[indexMCpiKK]->Fill(pt1);
671           fPtProng2Hist[indexMCpiKK]->Fill(pt2);
672         }
673       }
674      
675    
676     }
677     
678     Float_t tmp[26];
679     if(fFillNtuple>0){
680       
681       if ((fFillNtuple==1 && (isPhiKKpi || isPhipiKK)) || (fFillNtuple==2 && (isK0starKKpi || isK0starpiKK)) || (fFillNtuple==3 && (isKKpi || ispiKK))){
682         
683         tmp[0]=Float_t(labDs);
684         tmp[1]=Float_t(retCodeAnalysisCuts);
685         tmp[2]=Float_t(pdgCode0);  
686         tmp[3]=d->PtProng(0);
687         tmp[4]=d->PtProng(1);
688         tmp[5]=d->PtProng(2);
689         tmp[6]=d->Pt();
690         tmp[7]=d->CosPointingAngle();
691         tmp[8]=d->DecayLength();
692         tmp[9]=d->Xv();
693         tmp[10]=d->Yv();
694         tmp[11]=d->Zv();
695         tmp[12]=d->InvMassDsKKpi();
696         tmp[13]=d->InvMassDspiKK();
697         tmp[14]=d->GetSigmaVert();
698         tmp[15]=d->Getd0Prong(0);
699         tmp[16]=d->Getd0Prong(1);
700         tmp[17]=d->Getd0Prong(2);
701         tmp[18]=d->GetDCA();
702         tmp[19]=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
703         tmp[20]=d->InvMass2Prongs(0,1,321,321);
704         tmp[21]=d->InvMass2Prongs(1,2,321,321);
705         tmp[22]=d->CosPiDsLabFrameKKpi();          
706         tmp[23]=d->CosPiDsLabFramepiKK();       
707         tmp[24]=d->CosPiKPhiRFrameKKpi();          
708         tmp[25]=d->CosPiKPhiRFramepiKK();       
709         
710         
711         fNtupleDs->Fill(tmp);
712         PostData(4,fNtupleDs);
713       }  
714     }
715     
716     if(unsetvtx) d->UnsetOwnPrimaryVtx();
717     if(recVtx)fProdCuts->CleanOwnPrimaryVtx(d,aod,origownvtx);
718   }
719  
720   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
721   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
722
723   PostData(1,fOutput); 
724   PostData(3,fCounter);    
725
726   return;
727 }
728
729 //_________________________________________________________________
730
731 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
732 {
733   // Terminate analysis
734   //
735   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
736   fOutput = dynamic_cast<TList*> (GetOutputData(1));
737   if (!fOutput) {     
738     printf("ERROR: fOutput not available\n");
739     return;
740   }
741   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
742   if(fHistNEvents){
743     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
744   }else{
745     printf("ERROR: fHistNEvents not available\n");
746     return;
747   }
748   return;
749 }
750