]> 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 in histo binning
[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
c1a78ad6 278 fHistNtrCorrEvSel->Sumw2();
279 fHistNtrCorrEvWithCand->Sumw2();
280 fHistNtrCorrEvWithD->Sumw2();
0c514eec 281 fHistGenPrimaryParticlesInelGt0->Sumw2();
c1a78ad6 282 fOutput->Add(fHistNtrCorrEvSel);
283 fOutput->Add(fHistNtrCorrEvWithCand);
284 fOutput->Add(fHistNtrCorrEvWithD);
285 fOutput->Add(fHistNtrEta16vsNtrEta1);
286 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
287 fOutput->Add(fHistNtrVsZvtx);
288 fOutput->Add(fHistNtrCorrVsZvtx);
43093cfa 289
11457f09 290 fOutput->Add(fHistNtrVsNchMC);
291 fOutput->Add(fHistNtrCorrVsNchMC);
292 fOutput->Add(fHistNtrVsNchMCPrimary);
293 fOutput->Add(fHistNtrCorrVsNchMCPrimary);
294 fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
295 fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
0c514eec 296 fOutput->Add(fHistGenPrimaryParticlesInelGt0);
511c4e3e 297 fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
11457f09 298
43093cfa 299
300 fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
301 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
302 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
303 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
304 fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
305 fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
306 fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
307 fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
308 fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
309 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
310 fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
311 fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
312 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
313 fHistNEvents->Sumw2();
314 fHistNEvents->SetMinimum(0);
315 fOutput->Add(fHistNEvents);
316
f73a90da 317 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 318
f73a90da 319 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 320
7826c36d 321 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 322
f73a90da 323 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 324
f73a90da 325 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 326
327 fOutput->Add(fPtVsMassVsMult);
328 fOutput->Add(fPtVsMassVsMultUncorr);
329 fOutput->Add(fPtVsMassVsMultNoPid);
330 fOutput->Add(fPtVsMassVsMultPart);
331 fOutput->Add(fPtVsMassVsMultAntiPart);
332
333 if(fDoImpPar) CreateImpactParameterHistos();
334
335 fCounter = new AliNormalizationCounter("NormCounterCorrMult");
336 fCounter->SetStudyMultiplicity(kTRUE,1.);
337 fCounter->Init();
338
339 fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
340 fCounterU->SetStudyMultiplicity(kTRUE,1.);
341 fCounterU->Init();
342
343 fOutputCounters = new TList();
344 fOutputCounters->SetOwner();
345 fOutputCounters->SetName("OutputCounters");
346 fOutputCounters->Add(fCounter);
347 fOutputCounters->Add(fCounterU);
348
349 PostData(1,fOutput);
350 PostData(2,fListCuts);
351 PostData(3,fOutputCounters);
6c06d0fe 352 PostData(4,fListProfiles);
353
354
43093cfa 355
356 return;
357}
358
359
360
361//________________________________________________________________________
362void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
363{
364 // Execute analysis for current event:
365 // heavy flavor candidates association to MC truth
366
43093cfa 367 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
368
369 // AliAODTracklets* tracklets = aod->GetTracklets();
370 //Int_t ntracklets = tracklets->GetNumberOfTracklets();
371
372
373 TClonesArray *arrayCand = 0;
374 TString arrayName="";
375 UInt_t pdgDau[3];
376 Int_t nDau=0;
377 Int_t selbit=0;
378 if(fPdgMeson==411){
379 arrayName="Charm3Prong";
380 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
381 nDau=3;
382 selbit=AliRDHFCuts::kDplusCuts;
383 }else if(fPdgMeson==421){
384 arrayName="D0toKpi";
385 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
386 nDau=2;
387 selbit=AliRDHFCuts::kD0toKpiCuts;
388 }else if(fPdgMeson==413){
389 arrayName="Dstar";
390 pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211;
391 nDau=3;
392 selbit=AliRDHFCuts::kDstarCuts;
393 }
6c8b22ea 394
43093cfa 395 if(!aod && AODEvent() && IsStandardAOD()) {
396 // In case there is an AOD handler writing a standard AOD, use the AOD
397 // event in memory rather than the input (ESD) event.
398 aod = dynamic_cast<AliAODEvent*> (AODEvent());
399 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
400 // have to taken from the AOD event hold by the AliAODExtension
401 AliAODHandler* aodHandler = (AliAODHandler*)
402 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
403 if(aodHandler->GetExtensions()) {
404 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
405 AliAODEvent *aodFromExt = ext->GetAOD();
406 arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
407 }
408 } else if(aod) {
409 arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
410 }
43093cfa 411
412 if(!aod || !arrayCand) {
413 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
414 return;
415 }
416
417
418
419 // fix for temporary bug in ESDfilter
420 // the AODs with null vertex pointer didn't pass the PhysSel
421 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
422
423 Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
424 Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
425
426 fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1);
427 fHistNEvents->Fill(0); // count event
428
43093cfa 429 Double_t countTreta1corr=countTreta1;
f1c9cf1c 430 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
431 if(vtx1){
432 if(vtx1->GetNContributors()>0){
433 fHistNEvents->Fill(1);
434 TProfile* estimatorAvg = GetEstimatorHistogram(aod);
435 if(estimatorAvg){
436 countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
437 }
43093cfa 438 }
439 }
440
c1a78ad6 441 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countTreta1corr);
43093cfa 442
443 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
444
445 if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
446 if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
c1a78ad6 447 if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
448 if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
43093cfa 449
43093cfa 450
c1a78ad6 451 if(!isEvSel)return;
452 fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTr);
453 fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
454 if(vtx1){
455 fHistNtrVsZvtx->Fill(vtx1->GetZ(),countTreta1);
456 fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countTreta1corr);
457 }
458
43093cfa 459 TClonesArray *arrayMC=0;
460 AliAODMCHeader *mcHeader=0;
461
462 // load MC particles
463 if(fReadMC){
464
465 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
466 if(!arrayMC) {
467 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
468 return;
469 }
470 // load MC header
471 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
472 if(!mcHeader) {
473 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
474 return;
475 }
11457f09 476
477
478 Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
479 Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
480 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
0c514eec 481 if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
482 fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary);
483 }
11457f09 484 fHistNtrVsNchMC->Fill(nChargedMC,countTreta1);
485 fHistNtrCorrVsNchMC->Fill(nChargedMC,countTreta1corr);
486
487 fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1);
488 fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1corr);
489
490 fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1);
491 fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1corr);
492
511c4e3e 493 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary);
43093cfa 494 }
495
496 Int_t nCand = arrayCand->GetEntriesFast();
d2cdcb07 497 Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
498 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
499 Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
500 Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
43093cfa 501
502 for (Int_t iCand = 0; iCand < nCand; iCand++) {
503 AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
504 fHistNEvents->Fill(7);
505 if(fUseBit && !d->HasSelectionBit(selbit)){
506 fHistNEvents->Fill(8);
507 continue;
508 }
509
510 Double_t ptCand = d->Pt();
511 Double_t rapid=d->Y(fPdgMeson);
512 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
a8f9e77f 513 if(!isFidAcc) continue;
588dc385 514 Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
43093cfa 515 Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
516 if(passTopolCuts==0) continue;
517 nSelectedNoPID++;
518 fHistNEvents->Fill(9);
519 if(passAllCuts){
520 nSelectedPID++;
521 fHistNEvents->Fill(10);
522 }
523 Bool_t isPrimary=kTRUE;
524 Int_t labD=-1;
525 Double_t trueImpParXY=9999.;
526 Double_t impparXY=d->ImpParXY()*10000.;
527 Double_t dlen=0.1; //FIXME
528 Double_t mass[2];
529 if(fPdgMeson==411){
530 mass[0]=d->InvMass(nDau,pdgDau);
531 mass[1]=-1.;
d2cdcb07 532 if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
43093cfa 533 }else if(fPdgMeson==421){
534 UInt_t pdgdaughtersD0[2]={211,321};//pi,K
535 UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
536 mass[0]=d->InvMass(2,pdgdaughtersD0);
537 mass[1]=d->InvMass(2,pdgdaughtersD0bar);
d2cdcb07 538 if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
43093cfa 539 }else if(fPdgMeson==413){
540 // FIXME
541 AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
542 mass[0]=temp->DeltaInvMass();
543 mass[1]=-1.;
d2cdcb07 544 if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
43093cfa 545 }
546 for(Int_t iHyp=0; iHyp<2; iHyp++){
547 if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
548 Double_t invMass=mass[iHyp];
549 Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
550
551 if(fReadMC){
552
553 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
554 Bool_t fillHisto=fDoImpPar;
555 if(labD>=0){
556 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
557 Int_t code=partD->GetPdgCode();
558 if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
559 if(code<0 && iHyp==0) fillHisto=kFALSE;
560 if(code>0 && iHyp==1) fillHisto=kFALSE;
561 if(!isPrimary){
562 if(fPdgMeson==411){
563 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
564 }else if(fPdgMeson==421){
7826c36d 565 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
43093cfa 566 }else if(fPdgMeson==413){
567 trueImpParXY=0.; /// FIXME
568 }
569 Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
a8f9e77f 570 if(fillHisto && passAllCuts){
43093cfa 571 fHistMassPtImpPar[2]->Fill(arrayForSparse);
572 fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
573 }
574 }else{
a8f9e77f 575 if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
43093cfa 576 }
577 }else{
a8f9e77f 578 if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
43093cfa 579 }
580 if(fPdgMeson==421){
581 if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
582 if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
583 }
584 }
585
f3fddcd9 586 if(fPdgMeson==421){
587 if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
588 if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
589 }
43093cfa 590
a8f9e77f 591 fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
592 if(passAllCuts){
593 fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
594 fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
595 // Add separation between part antipart
596 if(fPdgMeson==411){
597 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
598 else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
599 }else if(fPdgMeson==421){
600 if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
601 if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
602 }else if(fPdgMeson==413){
603 // FIXME ADD Dstar!!!!!!!!
604 }
605
606 if(fDoImpPar){
607 fHistMassPtImpPar[0]->Fill(arrayForSparse);
43093cfa 608 }
a8f9e77f 609
43093cfa 610 }
611
612 }
613 }
614 fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
615 fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
c1a78ad6 616 fHistNtrCorrEvSel->Fill(countTreta1corr);
617 if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countTreta1corr);
618 if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countTreta1corr);
d2cdcb07 619
43093cfa 620 PostData(1,fOutput);
621 PostData(2,fListCuts);
f3fddcd9 622 PostData(3,fOutputCounters);
6c06d0fe 623
43093cfa 624 return;
625}
43093cfa 626//________________________________________________________________________
627void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
628 // Histos for impact paramter study
7826c36d 629 // mass . pt , impact parameter , decay length , multiplicity
43093cfa 630
7826c36d 631 Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
632 Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
633 Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
43093cfa 634
635 fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
636 "Mass vs. pt vs.imppar - All",
637 5,nbins,xmin,xmax);
638 fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
639 "Mass vs. pt vs.imppar - promptD",
640 5,nbins,xmin,xmax);
641 fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
642 "Mass vs. pt vs.imppar - DfromB",
643 5,nbins,xmin,xmax);
644 fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
645 "Mass vs. pt vs.true imppar -DfromB",
646 5,nbins,xmin,xmax);
647 fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
648 "Mass vs. pt vs.imppar - backgr.",
649 5,nbins,xmin,xmax);
650 for(Int_t i=0; i<5;i++){
651 fOutput->Add(fHistMassPtImpPar[i]);
652 }
653}
654
655//________________________________________________________________________
656void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
657{
658 // Terminate analysis
659 //
660 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
661
662 fOutput = dynamic_cast<TList*> (GetOutputData(1));
663 if (!fOutput) {
664 printf("ERROR: fOutput not available\n");
665 return;
666 }
2247bb0c 667
43093cfa 668 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
2247bb0c 669 if(!fHistNEvents){
670 printf("ERROR: fHistNEvents not available\n");
671 return;
672 }
c1a78ad6 673 printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
43093cfa 674
675 return;
676}
677//_________________________________________________________________________________________________
678Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
679 //
680 // checking whether the mother of the particles come from a charm or a bottom quark
681 //
682
683 Int_t pdgGranma = 0;
684 Int_t mother = 0;
685 mother = mcPartCandidate->GetMother();
686 Int_t istep = 0;
687 Int_t abspdgGranma =0;
688 Bool_t isFromB=kFALSE;
689 Bool_t isQuarkFound=kFALSE;
690 while (mother >0 ){
691 istep++;
692 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
693 if (mcGranma){
694 pdgGranma = mcGranma->GetPdgCode();
695 abspdgGranma = TMath::Abs(pdgGranma);
696 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
697 isFromB=kTRUE;
698 }
699 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
700 mother = mcGranma->GetMother();
701 }else{
702 AliError("Failed casting the mother particle!");
703 break;
704 }
705 }
706
707 if(isFromB) return 5;
708 else return 4;
709}
710
711
712
713//____________________________________________________________________________
714TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
715 // Get Estimator Histogram from period event->GetRunNumber();
716 //
717 // If you select SPD tracklets in |eta|<1 you should use type == 1
718 //
719
720 Int_t runNo = event->GetRunNumber();
721 Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
722 if(runNo>114930 && runNo<117223) period = 0;
723 if(runNo>119158 && runNo<120830) period = 1;
724 if(runNo>122373 && runNo<126438) period = 2;
725 if(runNo>127711 && runNo<130841) period = 3;
726 if(period<0 || period>3) return 0;
727
728 return fMultEstimatorAvg[period];
729}