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