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