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