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