1 /**************************************************************************
2 * Copyright(c) 1998-2009, 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 **************************************************************************/
15 //*************************************************************************
16 // Class AliAnalysisTaskSEBtoJPSItoEle
17 // AliAnalysisTaskSE for the reconstruction
18 // of heavy-flavour decay candidates
19 // Author: C.Di Giglio, carmelo.digiglio@ba.infn.it
20 //*************************************************************************
21 class AliAnalysisTaskSE;
28 #include <TClonesArray.h>
29 #include <TDatabasePDG.h>
31 #include "AliAnalysisManager.h"
32 #include "AliAODHandler.h"
33 #include "AliAODEvent.h"
34 #include "AliAODMCParticle.h"
35 #include "AliAODMCHeader.h"
36 #include "AliAODRecoDecayHF2Prong.h"
38 #include "AliAnalysisBtoJPSItoEle.h"
39 #include "AliAnalysisTaskSEBtoJPSItoEle.h"
41 ClassImp(AliAnalysisTaskSEBtoJPSItoEle)
43 //______________________________________________________________________________
44 AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle() :
49 fhDecayTimeMCjpsifromB(0),
55 fhCosThetaPointing(0),
61 // default constructor
63 //_________________________________________________________________________________________________
64 AliAnalysisTaskSEBtoJPSItoEle::AliAnalysisTaskSEBtoJPSItoEle(const char* name) :
65 AliAnalysisTaskSE(name),
69 fhDecayTimeMCjpsifromB(0),
75 fhCosThetaPointing(0),
81 // standard constructor
83 // Output slot #1 writes into a TList container
84 DefineOutput(1,TList::Class()); //My private output
86 //_________________________________________________________________________________________________
87 AliAnalysisTaskSEBtoJPSItoEle::~AliAnalysisTaskSEBtoJPSItoEle()
96 //_________________________________________________________________________________________________
97 void AliAnalysisTaskSEBtoJPSItoEle::Init()
101 if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::Init() \n");
106 //_________________________________________________________________________________________________
107 void AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects()
109 // Create the output container
111 if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle::UserCreateOutputObjects() \n");
113 // Several histograms are more conveniently managed in a TList
114 fOutput = new TList();
118 fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
119 fhDecayTimeMCjpsifromB->Sumw2();
120 fhDecayTimeMCjpsifromB->SetMinimum(0);
121 fOutput->Add(fhDecayTimeMCjpsifromB);
124 fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.);
125 fhDecayTime->Sumw2();
126 fhDecayTime->SetMinimum(0);
127 fOutput->Add(fhDecayTime);
129 fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; M [GeV]; Entries",100,0.,3.2);
131 fhInvMass->SetMinimum(0);
132 fOutput->Add(fhInvMass);
134 fhD0 = new TH1F("fhD0", "Impact parameter; D0 [#mu m]; Entries",100,-5000.,5000.);
139 fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.);
141 fhD0D0->SetMinimum(0);
142 fOutput->Add(fhD0D0);
144 fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2);
145 fhCosThetaStar->Sumw2();
146 fhCosThetaStar->SetMinimum(0);
147 fOutput->Add(fhCosThetaStar);
149 fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1);
150 fhCosThetaPointing->Sumw2();
151 fhCosThetaPointing->SetMinimum(0);
152 fOutput->Add(fhCosThetaPointing);
154 fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000);
155 fhCtsVsD0D0->Sumw2();
156 fhCtsVsD0D0->SetMinimum(0);
157 fOutput->Add(fhCtsVsD0D0);
159 fNtupleJPSI = new TNtuple("fNtupleJPSI","J/#Psi pseudo-proper decay time & invariant mass","Xdecaytime:Mass");
160 fOutput->Add(fNtupleJPSI);
164 //_________________________________________________________________________________________________
165 void AliAnalysisTaskSEBtoJPSItoEle::UserExec(Option_t */*option*/)
167 // Execute analysis for current event:
169 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
171 TClonesArray *inputArrayJPSItoEle = 0;
173 if(!aod && AODEvent() && IsStandardAOD()) {
174 // In case there is an AOD handler writing a standard AOD, use the AOD
175 // event in memory rather than the input (ESD) event.
176 aod = dynamic_cast<AliAODEvent*> (AODEvent());
177 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
178 // have to taken from the AOD event hold by the AliAODExtension
179 AliAODHandler* aodHandler = (AliAODHandler*)
180 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
181 if(aodHandler->GetExtensions()) {
182 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
183 AliAODEvent *aodFromExt = ext->GetAOD();
184 // load Jpsi candidates
185 inputArrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle");
188 // load Jpsi candidates
189 inputArrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
192 if(!inputArrayJPSItoEle) {
193 printf("AliAnalysisTaskSECompareHF::UserExec: JPSItoEle branch not found!\n");
197 // AOD primary vertex
198 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
201 TClonesArray* mcArray =
202 dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
203 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
204 AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
206 // Read AOD Monte Carlo information
207 if(fOkAODMC) ReadAODMCInfo(aod,inputArrayJPSItoEle);
209 // loop over J/Psi->ee candidates
210 Int_t nInJPSItoEle = inputArrayJPSItoEle->GetEntriesFast();
211 printf("Number of JPSI->e+e- candidates: %d\n",nInJPSItoEle);
213 for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
214 AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArrayJPSItoEle->UncheckedAt(iJPSItoEle);
216 Int_t mcLabel = d->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
219 Double_t vtxSec[3] = {0.,0.,0.};
220 Double_t vtxPrim[3] = {0.,0.,0.};
221 d->GetSecondaryVtx(vtxSec);
222 Bool_t unsetvtx=kFALSE;
223 if(!d->GetOwnPrimaryVtx()) {
224 d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
228 // here analyze only if cuts are passed
230 if(d->Pt() < fPtCuts[1] && d->Pt() > fPtCuts[0]){ // apply pt cut only
233 AliDebug(2,"No MC particle found");
237 fhInvMass->Fill(d->InvMassJPSIee());
238 fhD0->Fill(10000*d->ImpParXY());
239 fhD0D0->Fill(1e8*d->Prodd0d0());
240 fhCosThetaStar->Fill(d->CosThetaStarJPSI());
241 fhCtsVsD0D0->Fill(d->CosThetaStarJPSI(),1e8*d->Prodd0d0());
242 fhCosThetaPointing->Fill(d->CosPointingAngle());
244 // compute X variable
245 AliAODVertex* primVertex = d->GetOwnPrimaryVtx();
246 vtxPrim[0] = primVertex->GetX();
247 vtxPrim[1] = primVertex->GetY();
248 vtxPrim[2] = primVertex->GetZ();
249 Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(d->Px()) + (vtxSec[1]-vtxPrim[1])*(d->Py()))/d->Pt();
250 Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/d->Pt();
252 fhDecayTime->Fill(10000*pseudoProperDecayTime);
253 // Post the data already here
256 fNtupleJPSI->Fill(pseudoProperDecayTime,d->InvMassJPSIee()); // fill N-tuple with X,M values
260 } // end of JPSItoEle candidates selection according to cuts
262 if(unsetvtx) d->UnsetOwnPrimaryVtx();
264 }// end loop on JPSI to ele candidates
267 //_________________________________________________________________________________________________
268 void AliAnalysisTaskSEBtoJPSItoEle::Terminate(Option_t */*option*/)
271 // Terminate analysis
274 if(fDebug > 1) printf("AliAnalysisTaskSEBtoJPSItoEle: Terminate() \n");
276 fOutput = dynamic_cast<TList*> (GetOutputData(1));
278 printf("ERROR: fOutput not available\n");
282 if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
283 fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
284 fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
285 fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
286 fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
287 fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
288 fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
289 fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
290 fNtupleJPSI = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleJPSI"));
293 Double_t* pseudoProperValues=0x0;
294 Double_t* invMassValues=0x0;
296 fCdfFit = new AliAnalysisBtoJPSItoEle();
298 printf("+++\n+++ AliAnalysisBtoJPSItoEle object instantiated ---> OK\n+++\n");
299 fCdfFit->ReadCandidates(fNtupleJPSI,pseudoProperValues,invMassValues,ncand);
300 printf("+++\n+++ X and M vectors filled starting from N-tuple ---> OK\n+++\n");
301 printf("+++\n+++ Number of candidates ---> %d J/psi \n+++\n",ncand);
302 fCdfFit->SetPtBin(0);
303 printf("+++\n+++ Pt bin setted n. ---> %d \n+++\n",fCdfFit->GetPtBin());
305 if(fOkAODMC) {fCdfFit->CloneMCtemplate(fhDecayTimeMCjpsifromB);}
307 //SetPathMCfile(".");
308 TH1F* hLocal = OpenLocalMCFile(fNameMCfile,fCdfFit->GetPtBin());
309 fCdfFit->CloneMCtemplate(hLocal);
312 printf("+++\n+++ MC template histo copied ---> OK\n+++\n");
313 //fCdfFit->DoMinimization(pseudoProperValues,invMassValues,ncand);
314 fCdfFit->DoMinimization();
319 //_________________________________________________________________________________________________
320 void AliAnalysisTaskSEBtoJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
323 // apply Pt cuts (lower cut, upper cut)
325 for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
327 //_________________________________________________________________________________________________
328 void AliAnalysisTaskSEBtoJPSItoEle::SetCutsD0(const Double_t cuts[9])
331 // apply D0 and JPSI cuts
333 for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
335 //_________________________________________________________________________________________________
336 void AliAnalysisTaskSEBtoJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray)
339 // Reads MC information if needed (for example in case of parametrized PID)
343 TClonesArray* mcArray =
344 dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
345 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
346 AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
349 AliAODMCHeader* mcHeader =
350 (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
352 printf("AliAnalysisTaskSEBtoJPSItoEle::UserExec: MC AOD header branch not found!\n");
355 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
357 //Double_t rmax = 0.000005;
358 Double_t mcPrimVtx[3];
360 // get MC primary Vtx
361 mcHeader->GetVertex(mcPrimVtx);
363 Int_t nInCandidates = inputArray->GetEntriesFast();
364 printf("Number of Candidates for MC analysis: %d\n",nInCandidates);
366 Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
367 Int_t labJPSIdaugh0=0;
368 Int_t labJPSIdaugh1=0;
370 for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
372 AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
373 Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
375 Double_t vtxPrim[3] = {0.,0.,0.};
376 Double_t vtxSec[3] = {0.,0.,0.};
377 dd->GetSecondaryVtx(vtxSec);
378 Bool_t unsetvtx=kFALSE;
379 if(!dd->GetOwnPrimaryVtx()) {
380 dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
384 if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
386 AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
387 AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
388 lab0 = trk0->GetLabel();
389 lab1 = trk1->GetLabel();
391 AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
393 printf("no MC particle\n");
397 while(part0->GetMother()>=0) {
398 labMother=part0->GetMother();
399 part0 = (AliAODMCParticle*)mcArray->At(labMother);
401 printf("no MC mother particle\n");
404 pdgMother = TMath::Abs(part0->GetPdgCode());
405 if(pdgMother==443) {//this for JPSI
406 labJPSIdaugh0=labMother;
411 AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
413 printf("no MC particle\n");
417 while(part1->GetMother()>=0) {
418 labMother=part1->GetMother();
419 part1 = (AliAODMCParticle*)mcArray->At(labMother);
421 printf("no MC mother particle\n");
424 pdgMother = TMath::Abs(part1->GetPdgCode());
425 if(pdgMother==443) {//this for JPSI
426 labJPSIdaugh1=labMother;
431 if (mcLabelSec == -1)
433 AliDebug(2,"No MC particle found");
437 if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
438 AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
439 pdgJPSI = partJPSI->GetPdgCode();
440 Int_t pdaugh0 = partJPSI->GetDaughter(0);
441 Int_t pdaugh1 = partJPSI->GetDaughter(1);
442 AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
443 AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
444 pdg0 = partDaugh0->GetPdgCode();
445 pdg1 = partDaugh1->GetPdgCode();
446 if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
447 labJPSIMother = partJPSI->GetMother();
448 AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
449 pdgJPSIMother = partJPSIMother->GetPdgCode();
451 // chech if JPSI comes from B
452 if( pdgJPSIMother==511 || pdgJPSIMother==521 ||
453 pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
454 pdgJPSIMother==513 || pdgJPSIMother==523 ||
455 pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
456 pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
457 pdgJPSIMother==515 || pdgJPSIMother==525 ||
458 pdgJPSIMother==531 || pdgJPSIMother==10531 ||
459 pdgJPSIMother==533 || pdgJPSIMother==10533 ||
460 pdgJPSIMother==20533 || pdgJPSIMother==535 ||
461 pdgJPSIMother==541 || pdgJPSIMother==10541 ||
462 pdgJPSIMother==543 || pdgJPSIMother==10543 ||
463 pdgJPSIMother==20543 || pdgJPSIMother==545)
464 { //this is for MC JPSI -> ee from B-hadrons
466 fhInvMass->Fill(dd->InvMassJPSIee());
467 fhD0->Fill(10000*dd->ImpParXY());
468 fhD0D0->Fill(1e8*dd->Prodd0d0());
469 fhCosThetaStar->Fill(dd->CosThetaStarJPSI());
470 fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
471 fhCosThetaPointing->Fill(dd->CosPointingAngle());
473 // compute X variable
474 AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
475 vtxPrim[0] = primVertex->GetX();
476 vtxPrim[1] = primVertex->GetY();
477 vtxPrim[2] = primVertex->GetZ();
478 Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
479 Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
481 fhDecayTime->Fill(10000*pseudoProperDecayTime);
482 // Post the data already here
485 fNtupleJPSI->Fill(pseudoProperDecayTime,dd->InvMassJPSIee()); // fill N-tuple with X,M values
487 Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
488 Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
489 fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);
491 // Post the data already here
494 } //this is for MC JPSI -> ee from B-hadrons
495 } //this is for MC JPSI -> ee
497 } // end of check if JPSI is signal
502 //_________________________________________________________________________________________________
503 TH1F* AliAnalysisTaskSEBtoJPSItoEle::OpenLocalMCFile(TString datadir, Int_t binNum)
506 // Open a local file with MC x distribution for JPSI from B-hadron
509 TH1F* hCsiMC = new TH1F();
510 TString useFile = datadir.Data();
511 useFile.Append("/CsiMCfromKine_10PtBins.root");
512 TFile *f = new TFile(useFile);
514 char processCsiMCXHisto[1024];
515 if(binNum == 0) sprintf(processCsiMCXHisto,"hPseudoProper");
516 else sprintf(processCsiMCXHisto,"hPseudoProper%d",binNum);
517 hCsiMC = (TH1F*)f->Get(processCsiMCXHisto);
518 scale=1/hCsiMC->Integral();
519 hCsiMC->Scale(scale);
520 printf ("Opening Histo with Csi_MC(x) template with n. %f entries", hCsiMC->GetEntries());