1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //*************************************************************************
19 // Class AliAnalysisTaskSEDvsMultiplicity
20 // AliAnalysisTaskSE for the D meson vs. multiplcity analysis
21 // Authors: Renu Bala, Zaida Conesa del Valle, Francesco Prino
22 /////////////////////////////////////////////////////////////
24 #include <TClonesArray.h>
28 #include <TDatabasePDG.h>
32 #include <THnSparse.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"
50 ClassImp(AliAnalysisTaskSEDvsMultiplicity)
53 //________________________________________________________________________
54 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity():
61 fHistNtrEta16vsNtrEta1(0),
62 fHistNtrCorrEta1vsNtrRawEta1(0),
64 fHistNtrCorrVsZvtx(0),
66 fHistNtrCorrVsNchMC(0),
67 fHistNtrVsNchMCPrimary(0),
68 fHistNtrCorrVsNchMCPrimary(0),
69 fHistNtrVsNchMCPhysicalPrimary(0),
70 fHistNtrCorrVsNchMCPhysicalPrimary(0),
71 fHistGenPrimaryParticlesInelGt0(0),
72 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
73 fHistNtrUnCorrEvSel(0),
75 fHistNtrCorrEvWithCand(0),
76 fHistNtrCorrEvWithD(0),
78 fPtVsMassVsMultNoPid(0),
79 fPtVsMassVsMultUncorr(0),
80 fPtVsMassVsMultPart(0),
81 fPtVsMassVsMultAntiPart(0),
98 // Default constructor
99 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
100 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
103 //________________________________________________________________________
104 AliAnalysisTaskSEDvsMultiplicity::AliAnalysisTaskSEDvsMultiplicity(const char *name, Int_t pdgMeson,AliRDHFCuts *cuts):
105 AliAnalysisTaskSE(name),
111 fHistNtrEta16vsNtrEta1(0),
112 fHistNtrCorrEta1vsNtrRawEta1(0),
114 fHistNtrCorrVsZvtx(0),
116 fHistNtrCorrVsNchMC(0),
117 fHistNtrVsNchMCPrimary(0),
118 fHistNtrCorrVsNchMCPrimary(0),
119 fHistNtrVsNchMCPhysicalPrimary(0),
120 fHistNtrCorrVsNchMCPhysicalPrimary(0),
121 fHistGenPrimaryParticlesInelGt0(0),
122 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary(0),
123 fHistNtrUnCorrEvSel(0),
124 fHistNtrCorrEvSel(0),
125 fHistNtrCorrEvWithCand(0),
126 fHistNtrCorrEvWithD(0),
128 fPtVsMassVsMultNoPid(0),
129 fPtVsMassVsMultUncorr(0),
130 fPtVsMassVsMultPart(0),
131 fPtVsMassVsMultAntiPart(0),
133 fLowmasslimit(1.765),
135 fRDCutsAnalysis(cuts),
140 fLowerImpPar(-2000.),
141 fHigherImpPar(2000.),
149 // Standard constructor
151 for(Int_t i=0; i<5; i++) fHistMassPtImpPar[i]=0;
152 for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
154 fNMassBins=200; // FIXME
155 SetMassLimits(0.,0.2); // FIXME
158 SetMassLimits(fPdgMeson,0.1);
160 // Default constructor
161 // Output slot #1 writes into a TList container
162 DefineOutput(1,TList::Class()); //My private output
163 // Output slot #2 writes cut to private output
164 DefineOutput(2,TList::Class());
165 // Output slot #3 writes cut to private output
166 DefineOutput(3,TList::Class());
167 // Output slot #4 writes cut to private output
168 DefineOutput(4,TList::Class());
170 //________________________________________________________________________
171 AliAnalysisTaskSEDvsMultiplicity::~AliAnalysisTaskSEDvsMultiplicity()
179 delete fListProfiles;
180 delete fRDCutsAnalysis;
183 for(Int_t i=0; i<4; i++) delete fMultEstimatorAvg[i];
184 for(Int_t i=0; i<5; i++){
185 delete fHistMassPtImpPar[i];
189 //_________________________________________________________________
190 void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Double_t lowlimit, Double_t uplimit){
191 // set invariant mass limits
192 if(uplimit>lowlimit){
193 fLowmasslimit = lowlimit;
194 fUpmasslimit = uplimit;
196 AliError("Wrong mass limits: upper value should be larger than lower one");
199 //_________________________________________________________________
200 void AliAnalysisTaskSEDvsMultiplicity::SetMassLimits(Int_t pdg, Double_t range){
201 // set invariant mass limits
202 Double_t mass=TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
203 SetMassLimits(mass-range,mass+range);
205 //________________________________________________________________________
206 void AliAnalysisTaskSEDvsMultiplicity::Init(){
210 printf("AnalysisTaskSEDvsMultiplicity::Init() \n");
212 fListCuts=new TList();
215 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCutsAnalysis)));
216 copycut->SetName("AnalysisCutsDplus");
217 fListCuts->Add(copycut);
218 }else if(fPdgMeson==421){
219 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCutsAnalysis)));
220 copycut->SetName("AnalysisCutsDzero");
221 fListCuts->Add(copycut);
222 }else if(fPdgMeson==413){
223 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCutsAnalysis)));
224 copycut->SetName("AnalysisCutsDStar");
225 fListCuts->Add(copycut);
227 PostData(2,fListCuts);
229 fListProfiles = new TList();
230 fListProfiles->SetOwner();
231 TString period[4]={"LHC10b","LHC10c","LHC10d","LHC10e"};
232 for(Int_t i=0; i<4; i++){
233 if(fMultEstimatorAvg[i]){
234 TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
235 hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
236 fListProfiles->Add(hprof);
239 PostData(4,fListProfiles);
244 //________________________________________________________________________
245 void AliAnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects()
247 // Create the output container
249 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity::UserCreateOutputObjects() \n");
251 // Several histograms are more conveniently managed in a TList
252 fOutput = new TList();
254 fOutput->SetName("OutputHistos");
256 fHistNtrUnCorrEvSel = new TH1F("hNtrUnCorrEvSel","Uncorrected tracklets multiplicity for selected events; Tracklets ; Entries",200,-0.5,199.5);
257 fHistNtrCorrEvSel = new TH1F("hNtrCorrEvSel","Corrected tracklets multiplicity for selected events; Tracklets ; Entries",200,-0.5,199.5);
258 fHistNtrCorrEvWithCand = new TH1F("hNtrCorrEvWithCand", "Tracklets multiplicity for events with D candidates; Tracklets ; Entries",200,-0.5,199.5);// Total multiplicity
259 fHistNtrCorrEvWithD = new TH1F("hNtrCorrEvWithD", "Tracklets multiplicity for events with D in mass region ; Tracklets ; Entries",200,-0.5,199.5); //
260 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
261 fHistNtrCorrEta1vsNtrRawEta1 = new TH2F("hNtrCorrEta1vsNtrRawEta1","Corrected Eta1 vs Eta1.0; Ntracklets #eta<1.0 corrected; Ntracklets #eta<1",200,-0.5,199.5,200,-0.5,199.5); //eta 1.6 vs eta 1.0 histogram
262 fHistNtrVsZvtx = new TH2F("hNtrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,-0.5,199.5); //
263 fHistNtrCorrVsZvtx = new TH2F("hNtrCorrVsZvtx","Ntracklet vs VtxZ; VtxZ;N_{tracklet};",300,-15,15,200,-0.5,199.5); //
265 fHistNtrVsNchMC = new TH2F("hNtrVsNchMC","Ntracklet vs NchMC; Nch;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
266 fHistNtrCorrVsNchMC = new TH2F("hNtrCorrVsNchMC","Ntracklet vs Nch; Nch;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
268 fHistNtrVsNchMCPrimary = new TH2F("hNtrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch (Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
269 fHistNtrCorrVsNchMCPrimary = new TH2F("hNtrCorrVsNchMCPrimary","Ntracklet vs Nch (Primary); Nch(Primary) ;N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
271 fHistNtrVsNchMCPhysicalPrimary = new TH2F("hNtrVsNchMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
272 fHistNtrCorrVsNchMCPhysicalPrimary = new TH2F("hNtrCorrVsMCPhysicalPrimary","Ntracklet vs Nch (Physical Primary); Nch (Physical Primary);N_{tracklet};",200,-0.5,199.5,200,-0.5,199.5); //
274 fHistGenPrimaryParticlesInelGt0 = new TH1F("hGenPrimaryParticlesInelGt0","Multiplcity of generated charged particles ; Nparticles ; Entries",200,-0.5,199.5);
276 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary = new TH3F("fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary", "MC: Nch (Physical Primary) vs Nch (Primary) vs Nch (Generated); Nch (Generated); Nch (Primary); Nch (Physical Primary)",200,-0.5,199.5,200,-0.5,199.5,200,-0.5,199.5);
278 fHistNtrUnCorrEvSel->Sumw2();
279 fHistNtrCorrEvSel->Sumw2();
280 fHistNtrCorrEvWithCand->Sumw2();
281 fHistNtrCorrEvWithD->Sumw2();
282 fHistGenPrimaryParticlesInelGt0->Sumw2();
283 fOutput->Add(fHistNtrUnCorrEvSel);
284 fOutput->Add(fHistNtrCorrEvSel);
285 fOutput->Add(fHistNtrCorrEvWithCand);
286 fOutput->Add(fHistNtrCorrEvWithD);
287 fOutput->Add(fHistNtrEta16vsNtrEta1);
288 fOutput->Add(fHistNtrCorrEta1vsNtrRawEta1);
289 fOutput->Add(fHistNtrVsZvtx);
290 fOutput->Add(fHistNtrCorrVsZvtx);
292 fOutput->Add(fHistNtrVsNchMC);
293 fOutput->Add(fHistNtrCorrVsNchMC);
294 fOutput->Add(fHistNtrVsNchMCPrimary);
295 fOutput->Add(fHistNtrCorrVsNchMCPrimary);
296 fOutput->Add(fHistNtrVsNchMCPhysicalPrimary);
297 fOutput->Add(fHistNtrCorrVsNchMCPhysicalPrimary);
298 fOutput->Add(fHistGenPrimaryParticlesInelGt0);
299 fOutput->Add(fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary);
302 fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
303 fHistNEvents->GetXaxis()->SetBinLabel(1,"nEvents total");
304 fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with Z vertex");
305 fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents selected");
306 fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected due to trigger");
307 fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to phys sel");
308 fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vertex cuts");
309 fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to pileup");
310 fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
311 fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
312 fHistNEvents->GetXaxis()->SetBinLabel(10,"D after cuts (No PID)");
313 fHistNEvents->GetXaxis()->SetBinLabel(11,"D after cuts + PID)");
314 fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
315 fHistNEvents->Sumw2();
316 fHistNEvents->SetMinimum(0);
317 fOutput->Add(fHistNEvents);
319 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.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
321 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.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
323 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.);
325 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.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
327 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.5,199.5,fNMassBins,fLowmasslimit,fUpmasslimit,48,0.,24.);
329 fOutput->Add(fPtVsMassVsMult);
330 fOutput->Add(fPtVsMassVsMultUncorr);
331 fOutput->Add(fPtVsMassVsMultNoPid);
332 fOutput->Add(fPtVsMassVsMultPart);
333 fOutput->Add(fPtVsMassVsMultAntiPart);
335 if(fDoImpPar) CreateImpactParameterHistos();
337 fCounter = new AliNormalizationCounter("NormCounterCorrMult");
338 fCounter->SetStudyMultiplicity(kTRUE,1.);
341 fCounterU = new AliNormalizationCounter("NormCounterUnCorrMult");
342 fCounterU->SetStudyMultiplicity(kTRUE,1.);
345 fOutputCounters = new TList();
346 fOutputCounters->SetOwner();
347 fOutputCounters->SetName("OutputCounters");
348 fOutputCounters->Add(fCounter);
349 fOutputCounters->Add(fCounterU);
352 PostData(2,fListCuts);
353 PostData(3,fOutputCounters);
354 PostData(4,fListProfiles);
363 //________________________________________________________________________
364 void AliAnalysisTaskSEDvsMultiplicity::UserExec(Option_t */*option*/)
366 // Execute analysis for current event:
367 // heavy flavor candidates association to MC truth
369 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
371 // AliAODTracklets* tracklets = aod->GetTracklets();
372 //Int_t ntracklets = tracklets->GetNumberOfTracklets();
375 TClonesArray *arrayCand = 0;
376 TString arrayName="";
381 arrayName="Charm3Prong";
382 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=211;
384 selbit=AliRDHFCuts::kDplusCuts;
385 }else if(fPdgMeson==421){
387 pdgDau[0]=211; pdgDau[1]=321; pdgDau[2]=0;
389 selbit=AliRDHFCuts::kD0toKpiCuts;
390 }else if(fPdgMeson==413){
392 pdgDau[0]=321; pdgDau[1]=211; pdgDau[2]=211;
394 selbit=AliRDHFCuts::kDstarCuts;
397 if(!aod && AODEvent() && IsStandardAOD()) {
398 // In case there is an AOD handler writing a standard AOD, use the AOD
399 // event in memory rather than the input (ESD) event.
400 aod = dynamic_cast<AliAODEvent*> (AODEvent());
401 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
402 // have to taken from the AOD event hold by the AliAODExtension
403 AliAODHandler* aodHandler = (AliAODHandler*)
404 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
405 if(aodHandler->GetExtensions()) {
406 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
407 AliAODEvent *aodFromExt = ext->GetAOD();
408 arrayCand=(TClonesArray*)aodFromExt->GetList()->FindObject(arrayName.Data());
411 arrayCand=(TClonesArray*)aod->GetList()->FindObject(arrayName.Data());
414 if(!aod || !arrayCand) {
415 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: Charm3Prong branch not found!\n");
421 // fix for temporary bug in ESDfilter
422 // the AODs with null vertex pointer didn't pass the PhysSel
423 if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
425 Int_t countTreta1=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
426 Int_t countTr=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.6,1.6);
428 fCounterU->StoreEvent(aod,fRDCutsAnalysis,fReadMC,countTreta1);
429 fHistNEvents->Fill(0); // count event
431 Double_t countTreta1corr=countTreta1;
432 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
434 if(vtx1->GetNContributors()>0){
435 fHistNEvents->Fill(1);
436 TProfile* estimatorAvg = GetEstimatorHistogram(aod);
438 countTreta1corr=AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,countTreta1,vtx1->GetZ(),fRefMult);
443 fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC,(Int_t)countTreta1corr);
445 Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
447 if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(3);
448 if(fRDCutsAnalysis->GetWhyRejection()==7) fHistNEvents->Fill(4);
449 if(fRDCutsAnalysis->GetWhyRejection()==6) fHistNEvents->Fill(5);
450 if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(6);
454 fHistNtrEta16vsNtrEta1->Fill(countTreta1,countTr);
455 fHistNtrCorrEta1vsNtrRawEta1->Fill(countTreta1,countTreta1corr);
457 fHistNtrVsZvtx->Fill(vtx1->GetZ(),countTreta1);
458 fHistNtrCorrVsZvtx->Fill(vtx1->GetZ(),countTreta1corr);
461 TClonesArray *arrayMC=0;
462 AliAODMCHeader *mcHeader=0;
467 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
469 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC particles branch not found!\n");
473 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
475 printf("AliAnalysisTaskSEDvsMultiplicity::UserExec: MC header branch not found!\n");
480 Int_t nChargedMC=AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(arrayMC,-1.0,1.0);
481 Int_t nChargedMCPrimary=AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(arrayMC,-1.0,1.0);
482 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(arrayMC,-1.0,1.0);
483 if(nChargedMCPhysicalPrimary>0){ // INEL>0 for |eta|<1
484 fHistGenPrimaryParticlesInelGt0->Fill(nChargedMCPhysicalPrimary);
486 fHistNtrVsNchMC->Fill(nChargedMC,countTreta1);
487 fHistNtrCorrVsNchMC->Fill(nChargedMC,countTreta1corr);
489 fHistNtrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1);
490 fHistNtrCorrVsNchMCPrimary->Fill(nChargedMCPrimary,countTreta1corr);
492 fHistNtrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1);
493 fHistNtrCorrVsNchMCPhysicalPrimary->Fill(nChargedMCPhysicalPrimary,countTreta1corr);
495 fHistNchMCVsNchMCPrimaryVsNchMCPhysicalPrimary->Fill(nChargedMC,nChargedMCPrimary,nChargedMCPhysicalPrimary);
498 Int_t nCand = arrayCand->GetEntriesFast();
499 Int_t nSelectedNoPID=0,nSelectedPID=0,nSelectedInMassPeak=0;
500 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
501 Double_t mDplusPDG = TDatabasePDG::Instance()->GetParticle(411)->Mass();
502 Double_t mDstarPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
504 for (Int_t iCand = 0; iCand < nCand; iCand++) {
505 AliAODRecoDecayHF *d = (AliAODRecoDecayHF*)arrayCand->UncheckedAt(iCand);
506 fHistNEvents->Fill(7);
507 if(fUseBit && !d->HasSelectionBit(selbit)){
508 fHistNEvents->Fill(8);
512 Double_t ptCand = d->Pt();
513 Double_t rapid=d->Y(fPdgMeson);
514 Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
515 if(!isFidAcc) continue;
516 Int_t passAllCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
517 Int_t passTopolCuts=fRDCutsAnalysis->GetIsSelectedCuts();
518 if(passTopolCuts==0) continue;
520 fHistNEvents->Fill(9);
523 fHistNEvents->Fill(10);
525 Bool_t isPrimary=kTRUE;
527 Double_t trueImpParXY=9999.;
528 Double_t impparXY=d->ImpParXY()*10000.;
529 Double_t dlen=0.1; //FIXME
532 mass[0]=d->InvMass(nDau,pdgDau);
534 if(TMath::Abs(mass[0]-mDplusPDG)<0.02) nSelectedInMassPeak++; //20 MeV for now... FIXME
535 }else if(fPdgMeson==421){
536 UInt_t pdgdaughtersD0[2]={211,321};//pi,K
537 UInt_t pdgdaughtersD0bar[2]={321,211};//K,pi
538 mass[0]=d->InvMass(2,pdgdaughtersD0);
539 mass[1]=d->InvMass(2,pdgdaughtersD0bar);
540 if(TMath::Abs(mass[0]-mD0PDG)<0.02 || TMath::Abs(mass[1]-mD0PDG)<0.02 ) nSelectedInMassPeak++; //20 MeV for now... FIXME
541 }else if(fPdgMeson==413){
543 AliAODRecoCascadeHF* temp = (AliAODRecoCascadeHF*)d;
544 mass[0]=temp->DeltaInvMass();
546 if(TMath::Abs(mass[0]-(mDstarPDG-mD0PDG))<0.001) nSelectedInMassPeak++; //1 MeV for now... FIXME
548 for(Int_t iHyp=0; iHyp<2; iHyp++){
549 if(mass[iHyp]<0.) continue; // for D+ and D* we have 1 mass hypothesis
550 Double_t invMass=mass[iHyp];
551 Double_t arrayForSparse[5]={invMass,ptCand,impparXY,dlen,countTreta1corr};
555 labD = d->MatchToMC(fPdgMeson,arrayMC,nDau,(Int_t*)pdgDau);
556 Bool_t fillHisto=fDoImpPar;
558 AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
559 Int_t code=partD->GetPdgCode();
560 if(CheckOrigin(arrayMC,partD)==5) isPrimary=kFALSE;
561 if(code<0 && iHyp==0) fillHisto=kFALSE;
562 if(code>0 && iHyp==1) fillHisto=kFALSE;
565 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDplus(mcHeader,arrayMC,partD)*10000.;
566 }else if(fPdgMeson==421){
567 trueImpParXY=AliVertexingHFUtils::GetTrueImpactParameterDzero(mcHeader,arrayMC,partD)*10000.;
568 }else if(fPdgMeson==413){
569 trueImpParXY=0.; /// FIXME
571 Double_t arrayForSparseTrue[5]={invMass,ptCand,trueImpParXY,dlen,countTreta1corr};
572 if(fillHisto && passAllCuts){
573 fHistMassPtImpPar[2]->Fill(arrayForSparse);
574 fHistMassPtImpPar[3]->Fill(arrayForSparseTrue);
577 if(fillHisto && passAllCuts) fHistMassPtImpPar[1]->Fill(arrayForSparse);
580 if(fillHisto && passAllCuts)fHistMassPtImpPar[4]->Fill(arrayForSparse);
583 if(TMath::Abs(labD)==fPdgMeson && fMCOption==2) continue;
584 if(TMath::Abs(labD)!=fPdgMeson && fMCOption==1) continue;
589 if(iHyp==0 && !(passTopolCuts&1)) continue; // candidate not passing as D0
590 if(iHyp==1 && !(passTopolCuts&2)) continue; // candidate not passing as D0bar
593 fPtVsMassVsMultNoPid->Fill(countTreta1corr,invMass,ptCand);
595 fPtVsMassVsMult->Fill(countTreta1corr,invMass,ptCand);
596 fPtVsMassVsMultUncorr->Fill(countTreta1,invMass,ptCand);
597 // Add separation between part antipart
599 if(d->GetCharge()>0) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
600 else fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
601 }else if(fPdgMeson==421){
602 if(passTopolCuts&1) fPtVsMassVsMultPart->Fill(countTreta1corr,invMass,ptCand);
603 if(passTopolCuts&2) fPtVsMassVsMultAntiPart->Fill(countTreta1corr,invMass,ptCand);
604 }else if(fPdgMeson==413){
605 // FIXME ADD Dstar!!!!!!!!
609 fHistMassPtImpPar[0]->Fill(arrayForSparse);
616 fCounter->StoreCandidates(aod,nSelectedNoPID,kTRUE);
617 fCounter->StoreCandidates(aod,nSelectedPID,kFALSE);
618 fHistNtrUnCorrEvSel->Fill(countTreta1);
619 fHistNtrCorrEvSel->Fill(countTreta1corr);
620 if(nSelectedPID>0) fHistNtrCorrEvWithCand->Fill(countTreta1corr);
621 if(nSelectedInMassPeak>0) fHistNtrCorrEvWithD->Fill(countTreta1corr);
624 PostData(2,fListCuts);
625 PostData(3,fOutputCounters);
629 //________________________________________________________________________
630 void AliAnalysisTaskSEDvsMultiplicity::CreateImpactParameterHistos(){
631 // Histos for impact paramter study
632 // mass . pt , impact parameter , decay length , multiplicity
634 Int_t nbins[5]={fNMassBins,200,fNImpParBins,50,100};
635 Double_t xmin[5]={fLowmasslimit,0.,fLowerImpPar,0.,0.};
636 Double_t xmax[5]={fUpmasslimit,20.,fHigherImpPar,1.,100.};
638 fHistMassPtImpPar[0]=new THnSparseF("hMassPtImpParAll",
639 "Mass vs. pt vs.imppar - All",
641 fHistMassPtImpPar[1]=new THnSparseF("hMassPtImpParPrompt",
642 "Mass vs. pt vs.imppar - promptD",
644 fHistMassPtImpPar[2]=new THnSparseF("hMassPtImpParBfeed",
645 "Mass vs. pt vs.imppar - DfromB",
647 fHistMassPtImpPar[3]=new THnSparseF("hMassPtImpParTrueBfeed",
648 "Mass vs. pt vs.true imppar -DfromB",
650 fHistMassPtImpPar[4]=new THnSparseF("hMassPtImpParBkg",
651 "Mass vs. pt vs.imppar - backgr.",
653 for(Int_t i=0; i<5;i++){
654 fOutput->Add(fHistMassPtImpPar[i]);
658 //________________________________________________________________________
659 void AliAnalysisTaskSEDvsMultiplicity::Terminate(Option_t */*option*/)
661 // Terminate analysis
663 if(fDebug > 1) printf("AnalysisTaskSEDvsMultiplicity: Terminate() \n");
665 fOutput = dynamic_cast<TList*> (GetOutputData(1));
667 printf("ERROR: fOutput not available\n");
671 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
673 printf("ERROR: fHistNEvents not available\n");
676 printf("Number of Analyzed Events = %d\n",(Int_t)fHistNEvents->GetBinContent(3));
680 //_________________________________________________________________________________________________
681 Int_t AliAnalysisTaskSEDvsMultiplicity::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
683 // checking whether the mother of the particles come from a charm or a bottom quark
688 mother = mcPartCandidate->GetMother();
690 Int_t abspdgGranma =0;
691 Bool_t isFromB=kFALSE;
692 Bool_t isQuarkFound=kFALSE;
695 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
697 pdgGranma = mcGranma->GetPdgCode();
698 abspdgGranma = TMath::Abs(pdgGranma);
699 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
702 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
703 mother = mcGranma->GetMother();
705 AliError("Failed casting the mother particle!");
710 if(isFromB) return 5;
716 //____________________________________________________________________________
717 TProfile* AliAnalysisTaskSEDvsMultiplicity::GetEstimatorHistogram(const AliVEvent* event){
718 // Get Estimator Histogram from period event->GetRunNumber();
720 // If you select SPD tracklets in |eta|<1 you should use type == 1
723 Int_t runNo = event->GetRunNumber();
724 Int_t period = -1; // 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
725 if(runNo>114930 && runNo<117223) period = 0;
726 if(runNo>119158 && runNo<120830) period = 1;
727 if(runNo>122373 && runNo<126438) period = 2;
728 if(runNo>127711 && runNo<130841) period = 3;
729 if(period<0 || period>3) return 0;
731 return fMultEstimatorAvg[period];