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