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