]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix for par file creation on lion with clang (Dario Berzano)
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEDvsMultiplicity.cxx
CommitLineData
43093cfa 1/**************************************************************************
2 * Copyright(c) 1998-2008, 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// Class AliAnalysisTaskSEDvsMultiplicity
20// AliAnalysisTaskSE for the D meson vs. multiplcity analysis
21// Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
22/////////////////////////////////////////////////////////////
23
24#include <TClonesArray.h>
25#include <TCanvas.h>
26#include <TList.h>
27#include <TString.h>
28#include <TDatabasePDG.h>
29#include <TH1F.h>
30#include <TH2F.h>
31#include <TH3F.h>
32#include <THnSparse.h>
33#include <TProfile.h>
34#include "AliAnalysisManager.h"
35#include "AliRDHFCuts.h"
36#include "AliRDHFCutsDplustoKpipi.h"
37#include "AliRDHFCutsDStartoKpipi.h"
38#include "AliRDHFCutsD0toKpi.h"
39#include "AliAODHandler.h"
40#include "AliAODEvent.h"
41#include "AliAODVertex.h"
42#include "AliAODTrack.h"
43#include "AliAODRecoDecayHF.h"
44#include "AliAODRecoCascadeHF.h"
45#include "AliAnalysisVertexingHF.h"
46#include "AliAnalysisTaskSE.h"
47#include "AliAnalysisTaskSEDvsMultiplicity.h"
48#include "AliNormalizationCounter.h"
49#include "AliVertexingHFUtils.h"
50ClassImp(AliAnalysisTaskSEDvsMultiplicity)
51
52
53//________________________________________________________________________
54AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
55AliAnalysisTaskSE(),
56 fOutput(0),
57 fListCuts(0),
58 fOutputCounters(0),
6c06d0fe 59 fListProfiles(0),
43093cfa 60 fHistNEvents(0),
c1a78ad6 61 fHistNtrEta16vsNtrEta1(0),
62 fHistNtrCorrEta1vsNtrRawEta1(0),
63 fHistNtrVsZvtx(0),
64 fHistNtrCorrVsZvtx(0),
11457f09 65 fHistNtrVsNchMC(0),
66 fHistNtrCorrVsNchMC(0),
67 fHistNtrVsNchMCPrimary(0),
68 fHistNtrCorrVsNchMCPrimary(0),
69 fHistNtrVsNchMCPhysicalPrimary(0),
70 fHistNtrCorrVsNchMCPhysicalPrimary(0),
0c514eec 71 fHistGenPrimaryParticlesInelGt0(0),
511c4e3e 72 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
c1a78ad6 73 fHistNtrCorrEvSel(0),
74 fHistNtrCorrEvWithCand(0),
75 fHistNtrCorrEvWithD(0),
43093cfa 76 fPtVsMassVsMult(0),
77 fPtVsMassVsMultNoPid(0),
78 fPtVsMassVsMultUncorr(0),
79 fPtVsMassVsMultPart(0),
80 fPtVsMassVsMultAntiPart(0),
81 fUpmasslimit(1.965),
82 fLowmasslimit(1.765),
7826c36d 83 fNMassBins(200),
43093cfa 84 fRDCutsAnalysis(0),
85 fCounter(0),
86 fCounterU(0),
87 fDoImpPar(kFALSE),
88 fNImpParBins(400),
89 fLowerImpPar(-2000.),
90 fHigherImpPar(2000.),
91 fReadMC(kFALSE),
92 fMCOption(0),
93 fUseBit(kTRUE),
94 fRefMult(9.5),
95 fPdgMeson(411)
96{
97 // Default constructor
98 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
99 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
100}
101
102//________________________________________________________________________
7826c36d 103AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts):
43093cfa 104 AliAnalysisTaskSE(name),
105 fOutput(0),
106 fListCuts(0),
107 fOutputCounters(0),
6c06d0fe 108 fListProfiles(0),
43093cfa 109 fHistNEvents(0),
c1a78ad6 110 fHistNtrEta16vsNtrEta1(0),
111 fHistNtrCorrEta1vsNtrRawEta1(0),
112 fHistNtrVsZvtx(0),
113 fHistNtrCorrVsZvtx(0),
11457f09 114 fHistNtrVsNchMC(0),
115 fHistNtrCorrVsNchMC(0),
116 fHistNtrVsNchMCPrimary(0),
117 fHistNtrCorrVsNchMCPrimary(0),
118 fHistNtrVsNchMCPhysicalPrimary(0),
119 fHistNtrCorrVsNchMCPhysicalPrimary(0),
0c514eec 120 fHistGenPrimaryParticlesInelGt0(0),
511c4e3e 121 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
c1a78ad6 122 fHistNtrCorrEvSel(0),
123 fHistNtrCorrEvWithCand(0),
124 fHistNtrCorrEvWithD(0),
43093cfa 125 fPtVsMassVsMult(0),
126 fPtVsMassVsMultNoPid(0),
127 fPtVsMassVsMultUncorr(0),
128 fPtVsMassVsMultPart(0),
129 fPtVsMassVsMultAntiPart(0),
130 fUpmasslimit(1.965),
131 fLowmasslimit(1.765),
7826c36d 132 fNMassBins(200),
43093cfa 133 fRDCutsAnalysis(cuts),
134 fCounter(0),
135 fCounterU(0),
136 fDoImpPar(kFALSE),
137 fNImpParBins(400),
138 fLowerImpPar(-2000.),
139 fHigherImpPar(2000.),
140 fReadMC(kFALSE),
141 fMCOption(0),
142 fUseBit(kTRUE),
143 fRefMult(9.5),
7826c36d 144 fPdgMeson(pdgMeson)
43093cfa 145{
146 //
147 // Standard constructor
148 //
149 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
150 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
7826c36d 151 if(fPdgMeson==413){
152 fNMassBins=200; // FIXME
153 SetMassLimits(0.,0.2); // FIXME
154 }else{
155 fNMassBins=200;
156 SetMassLimits(fPdgMeson,0.1);
157 }
43093cfa 158 // Default constructor
159 // Output slot #1 writes into a TList container
160 DefineOutput(1,TList::Class()); //My private output
161 // Output slot #2 writes cut to private output
162 DefineOutput(2,TList::Class());
163 // Output slot #3 writes cut to private output
7826c36d 164 DefineOutput(3,TList::Class());
6c06d0fe 165 // Output slot #4 writes cut to private output
166 DefineOutput(4,TList::Class());
7826c36d 167}
43093cfa 168//________________________________________________________________________
169AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
170{
171 //
172 // Destructor
173 //
174 delete fOutput;
175 delete fHistNEvents;
176 delete fListCuts;
6c06d0fe 177 delete fListProfiles;
43093cfa 178 delete fRDCutsAnalysis;
179 delete fCounter;
180 delete fCounterU;
2d11aff8 181 for(Int_t i=0; i<4; i++) delete fMultEstimatorAvg[i];
43093cfa 182 for(Int_t i=0; i<5; i++){
183 delete fHistMassPtImpPar[i];
184 }
185}
186
187//_________________________________________________________________
7826c36d 188void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
43093cfa 189 // set invariant mass limits
190 if(uplimit>lowlimit){
a8f9e77f 191 fLowmasslimit = lowlimit;
192 fUpmasslimit = uplimit;
7826c36d 193 }else{
194 AliError("Wrong mass limits: upper value should be larger than lower one");
43093cfa 195 }
196}
7826c36d 197//_________________________________________________________________
198void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
199 // set invariant mass limits
200 Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
201 SetMassLimits(mass-range,mass+range);
43093cfa 202}
203//________________________________________________________________________
204void AliAnalysisTaskSEDvsMultiplicity::Init(){
205 //
206 // Initialization
207 //
208 printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
209
210 fListCuts=new TList();
211
212 if(fPdgMeson==411){
213 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
214 copycut->SetName("AnalysisCutsDplus");
215 fListCuts->Add(copycut);
216 }else if(fPdgMeson==421){
217 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
218 copycut->SetName("AnalysisCutsDzero");
219 fListCuts->Add(copycut);
220 }else if(fPdgMeson==413){
221 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
222 copycut->SetName("AnalysisCutsDStar");
223 fListCuts->Add(copycut);
224 }
225 PostData(2,fListCuts);
226
6c06d0fe 227 fListProfiles = new TList();
228 fListProfiles->SetOwner();
229 TString period[4]={"LHC10b","LHC10c","LHC10d","LHC10e"};
230 for(Int_t i=0; i<4; i++){
231 if(fMultEstimatorAvg[i]){
232 TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
233 hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
234 fListProfiles->Add(hprof);
235 }
236 }
237 PostData(4,fListProfiles);
238
43093cfa 239 return;
240}
241
242//________________________________________________________________________
243void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
244{
245 // Create the output container
246 //
247 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
248
249 // Several histograms are more conveniently managed in a TList
250 fOutput = new TList();
251 fOutput->SetOwner();
252 fOutput->SetName("OutputHistos");
253
c1a78ad6 254 fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Tracklets multiplicity for selected events; Tracklets ; Entries",200,0.,200.);
255 fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,0.,200.);// Total multiplicity
256 fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,0.,200.); //
257 fHistNtrEta16vsNtrEta1 = new TH2F("hNtrEta16vsNtrEta1","Uncorrected Eta1.6 vs Eta1.0; Ntracklets #eta<1.0; Ntracklets #eta<1.6",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram
258 fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",200,-0.,200.,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram
259 fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); //
260 fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,0,200.); //
11457f09 261
262 fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",200,0,200,200,0,200.); //
263 fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",200,0,200,200,0,200.); //
264
265 fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",200,0,200,200,0,200.); //
266 fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",200,0,200,200,0,200.); //
267
268 fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,0,200,200,0,200.); //
269 fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,0,200,200,0,200.); //
43093cfa 270
0c514eec 271 fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",200,-0.5,199.5);
43093cfa 272
511c4e3e 273 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",200,0.,200.,200,0.,200.,200,0.,200.);
274
c1a78ad6 275 fHistNtrCorrEvSel->Sumw2();
276 fHistNtrCorrEvWithCand->Sumw2();
277 fHistNtrCorrEvWithD->Sumw2();
0c514eec 278 fHistGenPrimaryParticlesInelGt0->Sumw2();
c1a78ad6 279 fOutput->Add(fHistNtrCorrEvSel);
280 fOutput->Add(fHistNtrCorrEvWithCand);
281 fOutput->Add(fHistNtrCorrEvWithD);
282 fOutput->Add(fHistNtrEta16vsNtrEta1);
283 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
284 fOutput->Add(fHistNtrVsZvtx);
285 fOutput->Add(fHistNtrCorrVsZvtx);
43093cfa 286
11457f09 287 fOutput->Add(fHistNtrVsNchMC);
288 fOutput->Add(fHistNtrCorrVsNchMC);
289 fOutput->Add(fHistNtrVsNchMCPrimary);
290 fOutput->Add(fHistNtrCorrVsNchMCPrimary);
291 fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
292 fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
0c514eec 293 fOutput->Add(fHistGenPrimaryParticlesInelGt0);
511c4e3e 294 fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
11457f09 295
43093cfa 296
297 fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
298 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
299 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
300 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
301 fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
302 fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
303 fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
304 fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
305 fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
306 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
307 fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
308 fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
309 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
310 fHistNEvents->Sumw2();
311 fHistNEvents->SetMinimum(0);
312 fOutput->Add(fHistNEvents);
313
7826c36d 314 fPtVsMassVsMult=new TH3F("hPtVsMassvsMult", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
43093cfa 315
7826c36d 316 fPtVsMassVsMultNoPid=new TH3F("hPtVsMassvsMultNoPid", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
43093cfa 317
7826c36d 318 fPtVsMassVsMultUncorr=new TH3F("hPtVsMassvsMultUncorr", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,-0.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
43093cfa 319
7826c36d 320 fPtVsMassVsMultPart=new TH3F("hPtVsMassvsMultPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
43093cfa 321
7826c36d 322 fPtVsMassVsMultAntiPart=new TH3F("hPtVsMassvsMultAntiPart", "D candidates: p_{t} vs mass vs tracklets multiplicity; Tracklets; Mass M [GeV/c^{2}]; p_{t} [GeV/c]",200,0.,200.,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
43093cfa 323
324 fOutput->Add(fPtVsMassVsMult);
325 fOutput->Add(fPtVsMassVsMultUncorr);
326 fOutput->Add(fPtVsMassVsMultNoPid);
327 fOutput->Add(fPtVsMassVsMultPart);
328 fOutput->Add(fPtVsMassVsMultAntiPart);
329
330 if(fDoImpPar) CreateImpactParameterHistos();
331
332 fCounter = new AliNormalizationCounter("NormCounterCorrMult");
333 fCounter->SetStudyMultiplicity(kTRUE,1.);
334 fCounter->Init();
335
336 fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
337 fCounterU->SetStudyMultiplicity(kTRUE,1.);
338 fCounterU->Init();
339
340 fOutputCounters = new TList();
341 fOutputCounters->SetOwner();
342 fOutputCounters->SetName("OutputCounters");
343 fOutputCounters->Add(fCounter);
344 fOutputCounters->Add(fCounterU);
345
346 PostData(1,fOutput);
347 PostData(2,fListCuts);
348 PostData(3,fOutputCounters);
6c06d0fe 349 PostData(4,fListProfiles);
350
351
43093cfa 352
353 return;
354}
355
356
357
358//________________________________________________________________________
359void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
360{
361 // Execute analysis for current event:
362 // heavy flavor candidates association to MC truth
363
43093cfa 364 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
365
366 // AliAODTracklets* tracklets = aod->GetTracklets();
367 //Int_t ntracklets = tracklets->GetNumberOfTracklets();
368
369
370 TClonesArray *arrayCand = 0;
371 TString arrayName="";
372 UInt_t pdgDau[3];
373 Int_t nDau=0;
374 Int_t selbit=0;
375 if(fPdgMeson==411){
376 arrayName="Charm3Prong";
377 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
378 nDau=3;
379 selbit=AliRDHFCuts::kDplusCuts;
380 }else if(fPdgMeson==421){
381 arrayName="D0toKpi";
382 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
383 nDau=2;
384 selbit=AliRDHFCuts::kD0toKpiCuts;
385 }else if(fPdgMeson==413){
386 arrayName="Dstar";
387 pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211;
388 nDau=3;
389 selbit=AliRDHFCuts::kDstarCuts;
390 }
6c8b22ea 391
43093cfa 392 if(!aod && AODEvent() && IsStandardAOD()) {
393 // In case there is an AOD handler writing a standard AOD, use the AOD
394 // event in memory rather than the input (ESD) event.
395 aod = dynamic_cast<AliAODEvent*> (AODEvent());
396 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
397 // have to taken from the AOD event hold by the AliAODExtension
398 AliAODHandler* aodHandler = (AliAODHandler*)
399 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
400 if(aodHandler->GetExtensions()) {
401 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
402 AliAODEvent *aodFromExt = ext->GetAOD();
403 arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
404 }
405 } else if(aod) {
406 arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
407 }
43093cfa 408
409 if(!aod || !arrayCand) {
410 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
411 return;
412 }
413
414
415
416 // fix for temporary bug in ESDfilter
417 // the AODs with null vertex pointer didn't pass the PhysSel
418 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
419
420 Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
421 Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
422
423 fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1);
424 fHistNEvents->Fill(0); // count event
425
43093cfa 426 Double_t countTreta1corr=countTreta1;
f1c9cf1c 427 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
428 if(vtx1){
429 if(vtx1->GetNContributors()>0){
430 fHistNEvents->Fill(1);
431 TProfile* estimatorAvg = GetEstimatorHistogram(aod);
432 if(estimatorAvg){
433 countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
434 }
43093cfa 435 }
436 }
437
c1a78ad6 438 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countTreta1corr);
43093cfa 439
440 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
441
442 if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
443 if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
c1a78ad6 444 if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
445 if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
43093cfa 446
43093cfa 447
c1a78ad6 448 if(!isEvSel)return;
449 fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTr);
450 fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
451 if(vtx1){
452 fHistNtrVsZvtx->Fill(vtx1->GetZ(),countTreta1);
453 fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countTreta1corr);
454 }
455
43093cfa 456 TClonesArray *arrayMC=0;
457 AliAODMCHeader *mcHeader=0;
458
459 // load MC particles
460 if(fReadMC){
461
462 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
463 if(!arrayMC) {
464 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
465 return;
466 }
467 // load MC header
468 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
469 if(!mcHeader) {
470 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
471 return;
472 }
11457f09 473
474
475 Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
476 Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
477 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
0c514eec 478 if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
479 fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary);
480 }
11457f09 481 fHistNtrVsNchMC->Fill(nChargedMC,countTreta1);
482 fHistNtrCorrVsNchMC->Fill(nChargedMC,countTreta1corr);
483
484 fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1);
485 fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1corr);
486
487 fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1);
488 fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1corr);
489
511c4e3e 490 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary);
43093cfa 491 }
492
493 Int_t nCand = arrayCand->GetEntriesFast();
d2cdcb07 494 Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
495 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
496 Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
497 Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
43093cfa 498
499 for (Int_t iCand = 0; iCand < nCand; iCand++) {
500 AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
501 fHistNEvents->Fill(7);
502 if(fUseBit && !d->HasSelectionBit(selbit)){
503 fHistNEvents->Fill(8);
504 continue;
505 }
506
507 Double_t ptCand = d->Pt();
508 Double_t rapid=d->Y(fPdgMeson);
509 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
a8f9e77f 510 if(!isFidAcc) continue;
588dc385 511 Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
43093cfa 512 Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
513 if(passTopolCuts==0) continue;
514 nSelectedNoPID++;
515 fHistNEvents->Fill(9);
516 if(passAllCuts){
517 nSelectedPID++;
518 fHistNEvents->Fill(10);
519 }
520 Bool_t isPrimary=kTRUE;
521 Int_t labD=-1;
522 Double_t trueImpParXY=9999.;
523 Double_t impparXY=d->ImpParXY()*10000.;
524 Double_t dlen=0.1; //FIXME
525 Double_t mass[2];
526 if(fPdgMeson==411){
527 mass[0]=d->InvMass(nDau,pdgDau);
528 mass[1]=-1.;
d2cdcb07 529 if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
43093cfa 530 }else if(fPdgMeson==421){
531 UInt_t pdgdaughtersD0[2]={211,321};//pi,K
532 UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
533 mass[0]=d->InvMass(2,pdgdaughtersD0);
534 mass[1]=d->InvMass(2,pdgdaughtersD0bar);
d2cdcb07 535 if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
43093cfa 536 }else if(fPdgMeson==413){
537 // FIXME
538 AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
539 mass[0]=temp->DeltaInvMass();
540 mass[1]=-1.;
d2cdcb07 541 if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
43093cfa 542 }
543 for(Int_t iHyp=0; iHyp<2; iHyp++){
544 if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
545 Double_t invMass=mass[iHyp];
546 Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
547
548 if(fReadMC){
549
550 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
551 Bool_t fillHisto=fDoImpPar;
552 if(labD>=0){
553 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
554 Int_t code=partD->GetPdgCode();
555 if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
556 if(code<0 && iHyp==0) fillHisto=kFALSE;
557 if(code>0 && iHyp==1) fillHisto=kFALSE;
558 if(!isPrimary){
559 if(fPdgMeson==411){
560 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
561 }else if(fPdgMeson==421){
7826c36d 562 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
43093cfa 563 }else if(fPdgMeson==413){
564 trueImpParXY=0.; /// FIXME
565 }
566 Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
a8f9e77f 567 if(fillHisto && passAllCuts){
43093cfa 568 fHistMassPtImpPar[2]->Fill(arrayForSparse);
569 fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
570 }
571 }else{
a8f9e77f 572 if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
43093cfa 573 }
574 }else{
a8f9e77f 575 if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
43093cfa 576 }
577 if(fPdgMeson==421){
578 if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
579 if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
580 }
581 }
582
f3fddcd9 583 if(fPdgMeson==421){
584 if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
585 if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
586 }
43093cfa 587
a8f9e77f 588 fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
589 if(passAllCuts){
590 fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
591 fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
592 // Add separation between part antipart
593 if(fPdgMeson==411){
594 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
595 else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
596 }else if(fPdgMeson==421){
597 if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
598 if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
599 }else if(fPdgMeson==413){
600 // FIXME ADD Dstar!!!!!!!!
601 }
602
603 if(fDoImpPar){
604 fHistMassPtImpPar[0]->Fill(arrayForSparse);
43093cfa 605 }
a8f9e77f 606
43093cfa 607 }
608
609 }
610 }
611 fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
612 fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
c1a78ad6 613 fHistNtrCorrEvSel->Fill(countTreta1corr);
614 if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countTreta1corr);
615 if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countTreta1corr);
d2cdcb07 616
43093cfa 617 PostData(1,fOutput);
618 PostData(2,fListCuts);
f3fddcd9 619 PostData(3,fOutputCounters);
6c06d0fe 620
43093cfa 621 return;
622}
43093cfa 623//________________________________________________________________________
624void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
625 // Histos for impact paramter study
7826c36d 626 // mass . pt , impact parameter , decay length , multiplicity
43093cfa 627
7826c36d 628 Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
629 Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
630 Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
43093cfa 631
632 fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
633 "Mass vs. pt vs.imppar - All",
634 5,nbins,xmin,xmax);
635 fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
636 "Mass vs. pt vs.imppar - promptD",
637 5,nbins,xmin,xmax);
638 fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
639 "Mass vs. pt vs.imppar - DfromB",
640 5,nbins,xmin,xmax);
641 fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
642 "Mass vs. pt vs.true imppar -DfromB",
643 5,nbins,xmin,xmax);
644 fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
645 "Mass vs. pt vs.imppar - backgr.",
646 5,nbins,xmin,xmax);
647 for(Int_t i=0; i<5;i++){
648 fOutput->Add(fHistMassPtImpPar[i]);
649 }
650}
651
652//________________________________________________________________________
653void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
654{
655 // Terminate analysis
656 //
657 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
658
659 fOutput = dynamic_cast<TList*> (GetOutputData(1));
660 if (!fOutput) {
661 printf("ERROR: fOutput not available\n");
662 return;
663 }
2247bb0c 664
43093cfa 665 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
2247bb0c 666 if(!fHistNEvents){
667 printf("ERROR: fHistNEvents not available\n");
668 return;
669 }
c1a78ad6 670 printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
43093cfa 671
672 return;
673}
674//_________________________________________________________________________________________________
675Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
676 //
677 // checking whether the mother of the particles come from a charm or a bottom quark
678 //
679
680 Int_t pdgGranma = 0;
681 Int_t mother = 0;
682 mother = mcPartCandidate->GetMother();
683 Int_t istep = 0;
684 Int_t abspdgGranma =0;
685 Bool_t isFromB=kFALSE;
686 Bool_t isQuarkFound=kFALSE;
687 while (mother >0 ){
688 istep++;
689 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
690 if (mcGranma){
691 pdgGranma = mcGranma->GetPdgCode();
692 abspdgGranma = TMath::Abs(pdgGranma);
693 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
694 isFromB=kTRUE;
695 }
696 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
697 mother = mcGranma->GetMother();
698 }else{
699 AliError("Failed casting the mother particle!");
700 break;
701 }
702 }
703
704 if(isFromB) return 5;
705 else return 4;
706}
707
708
709
710//____________________________________________________________________________
711TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
712 // Get Estimator Histogram from period event->GetRunNumber();
713 //
714 // If you select SPD tracklets in |eta|<1 you should use type == 1
715 //
716
717 Int_t runNo = event->GetRunNumber();
718 Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
719 if(runNo>114930 && runNo<117223) period = 0;
720 if(runNo>119158 && runNo<120830) period = 1;
721 if(runNo>122373 && runNo<126438) period = 2;
722 if(runNo>127711 && runNo<130841) period = 3;
723 if(period<0 || period>3) return 0;
724
725 return fMultEstimatorAvg[period];
726}