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