Updates from Chiara
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskSEDmesonsFilterCJ.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-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 //
17 // Task for selecting D mesons to be used as an input for D-Jet correlations
18 //
19 //-----------------------------------------------------------------------
20 // Authors:
21 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22 // A.Grelli (Utrecht University) a.grelli@uu.nl
23 // Xiaoming Zhang (LBNL)  XMZhang@lbl.gov
24 //-----------------------------------------------------------------------
25
26 #include <TDatabasePDG.h>
27 #include <TParticle.h>
28 #include <TVector3.h>
29 #include "TROOT.h"
30 #include <TH3F.h>
31
32 #include "AliRDHFCutsDStartoKpipi.h"
33 #include "AliRDHFCutsD0toKpi.h"
34 #include "AliAODMCHeader.h"
35 #include "AliAODHandler.h"
36 #include "AliAnalysisManager.h"
37 #include "AliLog.h"
38 #include "AliAODVertex.h"
39 #include "AliAODRecoDecay.h"
40 #include "AliAODRecoCascadeHF.h"
41 #include "AliAODRecoDecayHF2Prong.h"
42 #include "AliESDtrack.h"
43 #include "AliAODMCParticle.h"
44 #include "AliAnalysisTaskSEDmesonsFilterCJ.h"
45
46 ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
47
48 //_____________________________________________________________________________
49 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
50   AliAnalysisTaskSE(),
51   fUseMCInfo(kFALSE),
52   fCandidateType(0),
53   fCandidateName(""),
54   fPDGmother(0),
55   fNProngs(0),
56   fBranchName(""),
57   fOutput(0),
58 //fOutputCandidates(0),
59   fCuts(0),
60   fMinMass(0.),
61   fMaxMass(0.),
62   fCandidateArray(0)
63 //fIsSelectedArray(0)
64 {
65 //
66 // Default constructor
67 //
68
69   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
70   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
71 }
72
73 //_____________________________________________________________________________
74 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
75   AliAnalysisTaskSE(name),
76   fUseMCInfo(kFALSE),
77   fCandidateType(candtype),
78   fCandidateName(""),
79   fPDGmother(0),
80   fNProngs(0),
81 //fPDGdaughters(),
82   fBranchName(""),
83   fOutput(0),
84 //fOutputCandidates(0),
85   fCuts(cuts),
86   fMinMass(0.),
87   fMaxMass(0.),
88   fCandidateArray(0)
89 //fIsSelectedArray(0)
90 {
91 //
92 // Constructor. Initialization of Inputs and Outputs
93 //
94  
95   Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
96
97   for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
98   for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
99
100   const Int_t nptbins = fCuts->GetNPtBins();
101   Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };
102
103   switch (fCandidateType) {
104   case 0 :
105     fCandidateName = "D0";
106     fPDGmother = 421;
107     fNProngs = 2;
108     fPDGdaughters[0] = 211;  // pi 
109     fPDGdaughters[1] = 321;  // K
110     fPDGdaughters[2] = 0;    // empty
111     fPDGdaughters[3] = 0;    // empty
112     fBranchName = "D0toKpi";
113     break;
114   case 1 :
115     fCandidateName = "Dstar";
116     fPDGmother = 413;
117     fNProngs = 3;
118     fPDGdaughters[1] = 211; // pi soft
119     fPDGdaughters[0] = 421; // D0
120     fPDGdaughters[2] = 211; // pi fromD0
121     fPDGdaughters[3] = 321; // K from D0
122     fBranchName = "Dstar";
123
124     if (nptbins<=13) {
125       for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
126     } else {
127       AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
128     }
129     break;
130   default :
131     printf("%d not accepted!!\n",fCandidateType);
132     break;
133   }
134
135   if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
136   if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
137
138   DefineOutput(1, TList::Class());       // histos
139   DefineOutput(2, AliRDHFCuts::Class()); // my cuts
140 }
141
142 //_____________________________________________________________________________
143 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
144 {
145 //
146 // destructor
147 //
148
149   Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");  
150  
151   if (fOutput) { delete fOutput; fOutput = 0; }
152   if (fCuts)   { delete fCuts;   fCuts   = 0; }
153   if (fCandidateArray)  { delete fCandidateArray;  fCandidateArray  = 0; }
154 //if (fIsSelectedArray) { delete fIsSelectedArray; fIsSelectedArray = 0; }
155 }
156
157 //_____________________________________________________________________________
158 void AliAnalysisTaskSEDmesonsFilterCJ::Init()
159 {
160 //
161 // Initialization
162 //
163
164   if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
165
166   switch (fCandidateType) {
167     case 0: {
168               AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
169               copyfCutsDzero->SetName("AnalysisCutsDzero");
170               PostData(2, copyfCutsDzero);  // Post the data
171             } break;
172     case 1: {
173               AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
174               copyfCutsDstar->SetName("AnalysisCutsDStar");
175               PostData(2, copyfCutsDstar); // Post the cuts
176             } break;
177     default: return;
178   }
179
180   return;
181 }
182
183 //_____________________________________________________________________________
184 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
185
186 //
187 // output 
188 //
189
190   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
191   
192 //OpenFile(1);
193   fOutput = new TList(); fOutput->SetOwner();
194   DefineHistoForAnalysis(); // define histograms
195  
196 //fOutputCandidates= new TList();
197 //fOutputCandidates->SetOwner();
198
199   fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
200   fCandidateArray->SetName(Form("fCandidateArray%s",fCandidateName.Data()));
201
202 //fOutputCandidates->Add(fCandidateArray);
203 //fIsSelectedArray = new TClonesArray("Int_t",0);
204 //fIsSelectedArray->SetName(Form("fIsSelectedArray%s",suffix.Data()));
205 //fOutputCandidates->Add(fIsSelectedArray);
206
207   PostData(1, fOutput);
208 //PostData(3, fOutputCandidates);
209   return;
210 }
211
212 //_____________________________________________________________________________
213 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *)
214 {
215 //
216 // user exec
217 //
218
219   // add cadidate branch
220   fCandidateArray->Delete();
221   if (!(InputEvent()->FindListObject(Form("fCandidateArray%s",fCandidateName.Data())))) InputEvent()->AddObject(fCandidateArray);
222
223   // Load the event
224   AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
225  
226   TClonesArray *arrayDStartoD0pi = 0;
227   if (!aodEvent && AODEvent() && IsStandardAOD()) {
228
229     // In case there is an AOD handler writing a standard AOD, use the AOD 
230     // event in memory rather than the input (ESD) event.    
231     aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
232
233     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
234     // have to taken from the AOD event hold by the AliAODExtension
235     AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
236     if(aodHandler->GetExtensions()) {
237       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
238       AliAODEvent *aodFromExt = ext->GetAOD();
239       arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
240     }
241   } else {
242     arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
243   }
244
245   if (!arrayDStartoD0pi) {
246     AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
247     return;
248   } else {
249     AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   
250   }
251   
252   TClonesArray* mcArray = 0x0;
253   if (fUseMCInfo) {
254     mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
255     if (!mcArray) {
256       printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
257       return;
258     }
259   }
260
261   //Histograms
262   TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
263   TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
264
265   TH1F* hPtPion=0x0;
266   if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
267   hstat->Fill(0);
268
269   // fix for temporary bug in ESDfilter 
270   // the AODs with null vertex pointer didn't pass the PhysSel
271   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
272
273   //Event selection
274   Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
275 //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
276   if (!iseventselected) return;
277   hstat->Fill(1);
278
279   const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
280   hstat->Fill(2, nD);
281
282 //Int_t icountReco = 0;  // counters for efficiencies
283
284   //D* and D0 prongs needed to MatchToMC method
285   Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
286   Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
287
288   //D0 from D0 bar
289   Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K 
290   Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi 
291
292   Int_t iCand =0;
293   Int_t isSelected = 0;
294   AliAODRecoDecayHF *charmCand = 0;
295
296   Int_t mcLabel = -9999;
297   Int_t pdgMeson = 413;
298   if (fCandidateType==kD0toKpi) pdgMeson = 421;
299
300   for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
301     charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
302     if (!charmCand) return;
303
304
305     if (fUseMCInfo) { // Look in MC, try to simulate the z
306       if (fCandidateType==kDstartoKpipi) {
307         AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
308         mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
309       }
310
311       if (fCandidateType==kD0toKpi) mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
312
313       if (mcLabel<=0) continue;        
314     }
315
316     // region of interest + cuts
317     if (!fCuts->IsInFiducialAcceptance(charmCand->Pt(),charmCand->Y(pdgMeson))) continue;    
318
319     isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
320     if (!isSelected) continue;
321
322     new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
323 //  new ((*fIsSelectedArray)[iCand]) Int_t(isSelected);
324     iCand++;
325
326     Double_t masses[2];
327     Double_t ptD = charmCand->Pt();
328     if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
329
330       //softpion from D* decay
331       AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
332       AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();  
333       hstat->Fill(3);
334
335       // select D* in the D0 window.
336       // In the cut object window is loose to allow side bands
337       Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
338
339       // retrieve the corresponding pt bin for the candidate
340       // and set the expected D0 width (x3)
341 //    static const Int_t n = fCuts->GetNPtBins();
342       Int_t bin = fCuts->PtBin(ptD);
343       if(bin<0 || bin>=fCuts->GetNPtBins()) {
344         AliError(Form("Pt %.3f out of bounds",ptD));
345         continue;
346       }
347
348       AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
349       if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {   
350         masses[0] = temp->DeltaInvMass(); //D*
351         masses[1] = 0.; //dummy for D*
352
353         //D*  delta mass
354         hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
355
356         // D* pt and soft pion pt for good candidates           
357         Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
358         Double_t invmassDelta = temp->DeltaInvMass();
359         if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
360       }
361     } //Dstar specific
362
363     if (fCandidateType==kD0toKpi) { //D0->Kpi
364
365       //needed quantities
366       masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
367       masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
368       hstat->Fill(3);
369
370       // mass vs pt
371       if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
372       if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
373     } //D0 specific
374
375     charmCand = 0;  
376   } // end of D cand loop
377
378 //AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
379
380 //PostData(1,fOutput);
381 //PostData(3,fOutputCandidates);
382
383   return;
384 }
385
386 //_____________________________________________________________________________
387 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
388 {
389 // The Terminate() function is the last function to be called during
390 // a query. It always runs on the client, it can be used to present
391 // the results graphically or save the results to file.
392   
393   Info("Terminate"," terminate");
394   AliAnalysisTaskSE::Terminate();
395
396   fOutput = dynamic_cast<TList*>(GetOutputData(1));
397   if (!fOutput) {
398     printf("ERROR: fOutput not available\n");
399     return;
400   }
401
402   return;
403 }
404
405 //_____________________________________________________________________________
406 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
407 {
408 //
409 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
410 //
411
412   Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
413
414   // compute the Delta mass for the D*
415   if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
416 //cout << "mass ---------------" << endl;
417
418   fMinMass = mass - range;
419   fMaxMass = mass + range;
420
421   AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
422   if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
423
424   return;
425 }
426
427 //_____________________________________________________________________________
428 void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
429 {
430 //
431 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
432 //
433
434   if (uplimit>lowlimit) {
435     fMinMass = lowlimit;
436     fMaxMass = uplimit;
437   } else {
438     printf("Error! Lower limit larger than upper limit!\n");
439     fMinMass = uplimit - uplimit*0.2;
440     fMaxMass = uplimit;
441   }
442
443   return;
444 }
445
446 //_____________________________________________________________________________
447 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
448 {
449 //
450 // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
451 //
452
453   if (nptbins>30) {
454     AliInfo("Maximum number of bins allowed is 30!");
455     return kFALSE;
456   }
457
458   if (!width) return kFALSE;
459   for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
460
461   return kTRUE;
462 }
463
464
465 //_____________________________________________________________________________
466 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
467 {
468 //
469 // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
470 //
471
472   // Statistics 
473   TH1I* hstat = new TH1I("hstat","Statistics",4,-0.5,3.5);
474   hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
475   hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
476   hstat->GetXaxis()->SetBinLabel(3, "N cand");
477   hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
478 /*if(fUseMCInfo) {
479     hstat->GetXaxis()->SetBinLabel(7,"N D");
480     hstat->GetXaxis()->SetBinLabel(8,"N D in jet");
481   }*/
482   hstat->SetNdivisions(1);
483   fOutput->Add(hstat);
484
485   // Invariant mass related histograms
486   const Int_t nbinsmass = 200;
487   TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
488   hInvMass->SetStats(kTRUE);
489   hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
490   hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
491   fOutput->Add(hInvMass);
492   
493   if (fCandidateType==kDstartoKpipi) {
494     TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
495     hPtPion->SetStats(kTRUE);
496     hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
497     hPtPion->GetYaxis()->SetTitle("entries");
498     fOutput->Add(hPtPion);
499   }
500
501   return kTRUE; 
502 }