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