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