Dplus and Ds tasks use the new cuts classes (Francesco, Renu, Giacomo)
[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: AliITSCorrMapSDD.cxx 32906 2009-06-12 16:56:53Z prino $ */
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
46 ClassImp(AliAnalysisTaskSEDs)
47
48
49 //________________________________________________________________________
50 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs():
51   AliAnalysisTaskSE(),
52   fOutput(0), 
53   fHistNEvents(0),
54   fReadMC(kFALSE),
55   fNPtBins(0),
56   fListCuts(0),
57   fMassRange(0.2),
58   fProdCuts(0),
59   fAnalysisCuts(0)
60 {
61   // Default constructor
62 }
63
64 //________________________________________________________________________
65 AliAnalysisTaskSEDs::AliAnalysisTaskSEDs(const char *name, AliRDHFCutsDstoKKpi* productioncuts, AliRDHFCutsDstoKKpi* analysiscuts):
66   AliAnalysisTaskSE(name),
67   fOutput(0),
68   fHistNEvents(0),
69   fReadMC(kFALSE),
70   fNPtBins(0),
71   fListCuts(0),
72   fMassRange(0.2),
73   fProdCuts(productioncuts),
74   fAnalysisCuts(analysiscuts)
75 {
76   // Default constructor
77   // Output slot #1 writes into a TList container
78   Int_t nptbins=fAnalysisCuts->GetNPtBins();
79   Float_t *ptlim=fAnalysisCuts->GetPtBinLimits();
80   SetPtBins(nptbins,ptlim);
81  
82   DefineOutput(1,TList::Class());  //My private output
83
84   DefineOutput(2,TList::Class());
85 }
86
87 //________________________________________________________________________
88 void AliAnalysisTaskSEDs::SetPtBins(Int_t n, Float_t* lim){
89   // define pt bins for analysis
90   if(n>kMaxPtBins){
91     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
92     fNPtBins=kMaxPtBins;
93     fPtLimits[0]=0.;
94     fPtLimits[1]=1.;
95     fPtLimits[2]=3.;
96     fPtLimits[3]=5.;
97     fPtLimits[4]=10.;
98     for(Int_t i=5; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
99   }else{
100     fNPtBins=n;
101     for(Int_t i=0; i<fNPtBins+1; i++) fPtLimits[i]=lim[i];
102     for(Int_t i=fNPtBins+1; i<kMaxPtBins+1; i++) fPtLimits[i]=99999999.;
103   }
104   if(fDebug > 1){
105     printf("Number of Pt bins = %d\n",fNPtBins);
106     for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fPtLimits[i],fPtLimits[i+1]);    
107   }
108 }
109 //________________________________________________________________________
110 AliAnalysisTaskSEDs::~AliAnalysisTaskSEDs()
111 {
112   // Destructor
113   if (fOutput) {
114     delete fOutput;
115     fOutput = 0;
116   }
117   if (fListCuts) {
118     delete fListCuts;
119     fListCuts = 0;
120   }
121
122   if (fProdCuts) {
123     delete fProdCuts;
124     fProdCuts = 0;
125   }
126   if (fAnalysisCuts) {
127     delete fAnalysisCuts;
128     fAnalysisCuts = 0;
129   }
130 }  
131
132 //________________________________________________________________________
133 void AliAnalysisTaskSEDs::Init()
134 {
135   // Initialization
136
137   if(fDebug > 1) printf("AnalysisTaskSEDs::Init() \n");
138
139   fListCuts=new TList();
140   AliRDHFCutsDstoKKpi *production = new AliRDHFCutsDstoKKpi();
141   production=fProdCuts;
142   AliRDHFCutsDstoKKpi *analysis = new AliRDHFCutsDstoKKpi();
143   analysis=fAnalysisCuts;
144   
145   fListCuts->Add(production);
146   fListCuts->Add(analysis);
147   PostData(2,fListCuts);
148   return;
149 }
150
151 //________________________________________________________________________
152 void AliAnalysisTaskSEDs::UserCreateOutputObjects()
153 {
154   // Create the output container
155   //
156   if(fDebug > 1) printf("AnalysisTaskSEDs::UserCreateOutputObjects() \n");
157
158   // Several histograms are more conveniently managed in a TList
159   fOutput = new TList();
160   fOutput->SetOwner();
161   fOutput->SetName("OutputHistos");
162
163   Double_t massDs=TDatabasePDG::Instance()->GetParticle(431)->Mass();
164   Double_t minMass=massDs-0.5*fMassRange;
165   Double_t maxMass=massDs+0.5*fMassRange;
166   TString hisname;
167   Int_t index;
168   for(Int_t i=0;i<fNPtBins;i++){
169     index=GetHistoIndex(i);
170     hisname.Form("hMassAllPt%d",i);
171     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
172     fMassHist[index]->Sumw2();
173     hisname.Form("hMassAllPt%dCuts",i);
174     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
175     fMassHistCuts[index]->Sumw2();
176     hisname.Form("hCosPAllPt%d",i);
177     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
178     fCosPHist[index]->Sumw2();
179     hisname.Form("hDLenAllPt%d",i);
180     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
181     fDLenHist[index]->Sumw2();
182     hisname.Form("hDalitzAllPt%d",i);
183     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
184     fDalitz[index]->Sumw2();
185
186     index=GetSignalHistoIndex(i);    
187     hisname.Form("hMassSigPt%d",i);
188     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
189     fMassHist[index]->Sumw2();
190     hisname.Form("hMassSigPt%dCuts",i);
191     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
192     fMassHistCuts[index]->Sumw2();
193     hisname.Form("hCosPSigPt%d",i);
194     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
195     fCosPHist[index]->Sumw2();
196     hisname.Form("hDLenSigPt%d",i);
197     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
198     fDLenHist[index]->Sumw2();
199     hisname.Form("hDalitzSigPt%d",i);
200     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
201     fDalitz[index]->Sumw2();
202
203     hisname.Form("hMassBkgPt%d",i);
204     index=GetBackgroundHistoIndex(i);    
205     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
206     fMassHist[index]->Sumw2();
207     hisname.Form("hMassBkgPt%dCuts",i);
208     fMassHistCuts[index]=new TH1F(hisname.Data(),hisname.Data(),100,minMass,maxMass);
209     fMassHistCuts[index]->Sumw2();
210     hisname.Form("hCosPBkgPt%d",i);
211     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
212     fCosPHist[index]->Sumw2();
213     hisname.Form("hDLenBkgPt%d",i);
214     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
215     fDLenHist[index]->Sumw2();
216     hisname.Form("hDalitzBkgPt%d",i);
217     fDalitz[index]=new TH2F(hisname.Data(),hisname.Data(),100,0.,2.,100,0.,2.);
218     fDalitz[index]->Sumw2();
219   }
220
221   for(Int_t i=0; i<3*fNPtBins; i++){
222     fOutput->Add(fMassHist[i]);
223     fOutput->Add(fMassHistCuts[i]);
224     fOutput->Add(fCosPHist[i]);
225     fOutput->Add(fDLenHist[i]);
226     fOutput->Add(fDalitz[i]);
227   }
228
229   fChanHist[0] = new TH1F("hChanAll", "KKpi and piKK candidates",4,-0.5,3.5);
230   fChanHist[1] = new TH1F("hChanSig", "KKpi and piKK candidates",4,-0.5,3.5);
231   fChanHist[2] = new TH1F("hChanBkg", "KKpi and piKK candidates",4,-0.5,3.5);
232   fChanHistCuts[0] = new TH1F("hChanAllCuts", "KKpi and piKK candidates",4,-0.5,3.5);
233   fChanHistCuts[1] = new TH1F("hChanSigCuts", "KKpi and piKK candidates",4,-0.5,3.5);
234   fChanHistCuts[2] = new TH1F("hChanBkgCuts", "KKpi and piKK candidates",4,-0.5,3.5);
235   for(Int_t i=0;i<3;i++){
236     fChanHist[i]->Sumw2();
237     fChanHist[i]->SetMinimum(0);
238     fChanHistCuts[i]->Sumw2();
239     fChanHistCuts[i]->SetMinimum(0);
240     fOutput->Add(fChanHist[i]);
241     fOutput->Add(fChanHistCuts[i]);
242   }
243
244   fHistNEvents = new TH1F("hNEvents", "Number of processed events",3,-1.5,1.5);
245   fHistNEvents->Sumw2();
246   fHistNEvents->SetMinimum(0);
247   fOutput->Add(fHistNEvents);
248
249
250   return;
251 }
252
253 //________________________________________________________________________
254 void AliAnalysisTaskSEDs::UserExec(Option_t */*option*/)
255 {
256   // Ds selection for current event, fill mass histos and selecetion variable histo
257   // separate signal and backgound if fReadMC is activated
258
259   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
260   fHistNEvents->Fill(0); // count event
261   // Post the data already here
262   PostData(1,fOutput);
263   
264
265   TClonesArray *array3Prong = 0;
266   if(!aod && AODEvent() && IsStandardAOD()) {
267     // In case there is an AOD handler writing a standard AOD, use the AOD 
268     // event in memory rather than the input (ESD) event.    
269     aod = dynamic_cast<AliAODEvent*> (AODEvent());
270     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
271     // have to taken from the AOD event hold by the AliAODExtension
272     AliAODHandler* aodHandler = (AliAODHandler*) 
273       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
274     if(aodHandler->GetExtensions()) {
275       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
276       AliAODEvent *aodFromExt = ext->GetAOD();
277       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
278     }
279   } else {
280     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
281   }
282
283   if(!array3Prong) {
284     printf("AliAnalysisTaskSEDs::UserExec: Charm3Prong branch not found!\n");
285     return;
286   }
287
288  
289   TClonesArray *arrayMC=0;
290   AliAODMCHeader *mcHeader=0;
291
292   // AOD primary vertex
293   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
294   //    vtx1->Print();
295   
296   // load MC particles
297   if(fReadMC){
298     
299     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
300     if(!arrayMC) {
301       printf("AliAnalysisTaskSEDs::UserExec: MC particles branch not found!\n");
302       return;
303     }
304     
305   // load MC header
306     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
307     if(!mcHeader) {
308       printf("AliAnalysisTaskSEDs::UserExec: MC header branch not found!\n");
309       return;
310     }
311   }
312   
313   Int_t n3Prong = array3Prong->GetEntriesFast();
314   if(fDebug>1) printf("Number of Ds->KKpi: %d\n",n3Prong);
315   
316   
317   Int_t pdgDstoKKpi[3]={321,321,211};
318   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
319     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
320     
321     
322     Bool_t unsetvtx=kFALSE;
323     if(!d->GetOwnPrimaryVtx()){
324       d->SetOwnPrimaryVtx(vtx1);
325       unsetvtx=kTRUE;
326     }
327
328     Int_t retCodeProductionCuts=fProdCuts->IsSelected(d,AliRDHFCuts::kCandidate);
329     if(retCodeProductionCuts>0){
330       Int_t isKKpi=retCodeProductionCuts&1;
331       Int_t ispiKK=retCodeProductionCuts&2;
332 //       Int_t isPhi=retCodeProductionCuts&4;
333 //       Int_t isK0star=retCodeProductionCuts&8;
334       Double_t ptCand = d->Pt();
335       Int_t iPtBin=TMath::BinarySearch(fNPtBins,fPtLimits,(Float_t)ptCand);
336       Int_t retCodeAnalysisCuts=fAnalysisCuts->IsSelected(d,AliRDHFCuts::kCandidate);
337       Int_t isKKpiAC=retCodeAnalysisCuts&1;
338       Int_t ispiKKAC=retCodeAnalysisCuts&2;
339 //       Int_t isPhiAC=retCodeAnalysisCuts&4;
340 //       Int_t isK0starAC=retCodeAnalysisCuts&8;
341
342       Int_t labDs=-1;
343       if(fReadMC){
344         labDs = d->MatchToMC(431,arrayMC,3,pdgDstoKKpi);
345       }
346
347       Double_t dlen=d->DecayLength();
348       Double_t cosp=d->CosPointingAngle();
349       Int_t index=GetHistoIndex(iPtBin);
350       Int_t type=0;
351       if(isKKpi) type+=1;
352       if(ispiKK) type+=2;
353       Int_t typeAC=0;
354       if(isKKpiAC) typeAC+=1;
355       if(ispiKKAC) typeAC+=2;
356       fCosPHist[index]->Fill(cosp);
357       fDLenHist[index]->Fill(dlen);
358       fChanHist[0]->Fill(type);
359       if(retCodeAnalysisCuts>0) fChanHistCuts[0]->Fill(typeAC);
360       if(fReadMC){
361         if(labDs>=0) {    
362           index=GetSignalHistoIndex(iPtBin);
363           fCosPHist[index]->Fill(cosp);
364           fDLenHist[index]->Fill(dlen);
365           fChanHist[1]->Fill(type);       
366           if(retCodeAnalysisCuts>0) fChanHistCuts[1]->Fill(typeAC);
367         }else{
368           index=GetBackgroundHistoIndex(iPtBin);
369           fCosPHist[index]->Fill(cosp);
370           fDLenHist[index]->Fill(dlen);
371           fChanHist[2]->Fill(type);       
372           if(retCodeAnalysisCuts>0) fChanHistCuts[2]->Fill(typeAC);
373         }
374       }
375       if(isKKpi){
376         index=GetHistoIndex(iPtBin);
377         Double_t invMass=d->InvMassDsKKpi();
378         fMassHist[index]->Fill(invMass);
379         Double_t mass01=d->InvMass2Prongs(0,1,321,321);
380         Double_t mass12=d->InvMass2Prongs(1,2,321,211);
381         fDalitz[index]->Fill(mass01,mass12);
382         if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
383         if(fReadMC){
384           if(labDs>=0) {          
385             index=GetSignalHistoIndex(iPtBin);
386             fMassHist[index]->Fill(invMass);
387             fDalitz[index]->Fill(mass01,mass12);
388             if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
389           }else{
390             index=GetBackgroundHistoIndex(iPtBin);
391             fMassHist[index]->Fill(invMass);
392             fDalitz[index]->Fill(mass01,mass12);
393             if(retCodeAnalysisCuts>0 && isKKpiAC) fMassHistCuts[index]->Fill(invMass);
394           }
395         }       
396       }
397       if(ispiKK){
398         index=GetHistoIndex(iPtBin);
399         Double_t invMass=d->InvMassDspiKK();
400         fMassHist[index]->Fill(invMass);
401         if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
402         if(fReadMC){
403           if(labDs>=0) {          
404             index=GetSignalHistoIndex(iPtBin);
405             fMassHist[index]->Fill(invMass);
406             if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
407           }else{
408             index=GetBackgroundHistoIndex(iPtBin);
409             fMassHist[index]->Fill(invMass);
410             if(retCodeAnalysisCuts>0 && ispiKKAC) fMassHistCuts[index]->Fill(invMass);
411           }
412         }
413       }
414     }
415     if(unsetvtx) d->UnsetOwnPrimaryVtx();
416   }
417  
418    
419   PostData(1,fOutput);    
420   return;
421 }
422
423 //_________________________________________________________________
424
425 void AliAnalysisTaskSEDs::Terminate(Option_t */*option*/)
426 {
427   // Terminate analysis
428   //
429   if(fDebug > 1) printf("AnalysisTaskSEDs: Terminate() \n");
430   fOutput = dynamic_cast<TList*> (GetOutputData(1));
431   if (!fOutput) {     
432     printf("ERROR: fOutput not available\n");
433     return;
434   }
435   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
436   fChanHist[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAll"));
437   fChanHist[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSig"));
438   fChanHist[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkg"));
439   fChanHistCuts[0] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanAllCuts"));
440   fChanHistCuts[1] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanSigCuts"));
441   fChanHistCuts[2] = dynamic_cast<TH1F*>(fOutput->FindObject("hChanBkgCuts"));
442
443
444   TString hisname;
445   Int_t index;
446   for(Int_t i=0;i<fNPtBins;i++){
447
448     index=GetHistoIndex(i);
449     hisname.Form("hMassAllPt%d",i);
450     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
451     hisname.Form("hMassAllPt%dCuts",i);
452     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
453     hisname.Form("hCosPAllPt%d",i);
454     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
455     hisname.Form("hDLenAllPt%d",i);
456     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
457     hisname.Form("hDalitzAllPt%d",i);
458     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
459
460     index=GetSignalHistoIndex(i);    
461     hisname.Form("hMassSigPt%d",i);
462     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
463     hisname.Form("hMassSigPt%dCuts",i);
464     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
465     hisname.Form("hCosPSigPt%d",i);
466     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
467     hisname.Form("hDLenSigPt%d",i);
468     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
469     hisname.Form("hDalitzSigPt%d",i);
470     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
471
472     index=GetBackgroundHistoIndex(i);    
473     hisname.Form("hMassBkgPt%d",i);
474     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
475     hisname.Form("hMassBkgPt%dCuts",i);
476     fMassHistCuts[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
477     hisname.Form("hCosPBkgPt%d",i);
478     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
479     hisname.Form("hDLenBkgPt%d",i);
480     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
481     hisname.Form("hDalitzBkgPt%d",i);
482     fDalitz[index]=dynamic_cast<TH2F*>(fOutput->FindObject(hisname.Data()));
483
484   }
485   return;
486 }