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 AliAnalysisTaskSEJPSItoEle
17 // AliAnalysisTaskSE class to read both J/psi -> e+e- candidates and
18 // like-sign pair candidates and to store them into a specific
19 // stand-alone AOD file for J/psi into dieletrons analysis.
22 // 1) produces several histograms (including invariant mass distributions)
23 // for both unlike sign (US) and like sign (LS) pairs.
25 // 2) selects only J/Psi to dielectrons candidates and copies it to a
26 // specific AOD file, namely AliAODjpsi.root. The references
27 // to AliAODTrack objects are copied as well.
28 // The specific AliAODjpsi.root is thought as the input file
29 // to the AliAnalysisBtoJPSItoEle.h class, which performs (locally)
30 // the analysis of J/Psi from beauty hadrons decays.
32 // Author: C.Di Giglio, carmelo.digiglio@ba.infn.it
33 //*************************************************************************
34 class AliAnalysisTaskSE;
41 #include "TClonesArray.h"
42 #include "TDatabasePDG.h"
45 #include "TBits.h" //**********************
47 #include "AliAODEvent.h"
48 #include "AliAODMCParticle.h"
49 #include "AliAODMCHeader.h"
50 #include "AliAODRecoDecayHF2Prong.h"
51 #include "AliAODHandler.h"
52 #include "AliAnalysisManager.h"
53 #include "AliAnalysisTaskSEJPSItoEle.h"
54 #include "AliAODTrack.h"
55 #include "AliExternalTrackParam.h"
57 ClassImp(AliAnalysisTaskSEJPSItoEle)
59 //______________________________________________________________________________
60 AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle() :
63 fhDecayTimeMCjpsifromB(0),
71 fhCosThetaPointing(0),
94 // default constructor
96 //_________________________________________________________________________________________________
97 AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle(const char* name) :
98 AliAnalysisTaskSE(name),
100 fhDecayTimeMCjpsifromB(0),
108 fhCosThetaPointing(0),
121 fLsNormalization(1.),
124 fVerticesHFTClArr(0),
131 // standard constructor
133 // Input slot #0 works with an Ntuple
134 DefineInput(0, TChain::Class());
136 // Output slot #0 writes into a TTree container
137 // Output slot #1 writes into a TList container
138 DefineOutput(0, TTree::Class());
139 DefineOutput(1,TList::Class()); //My private output
141 //_________________________________________________________________________________________________
142 AliAnalysisTaskSEJPSItoEle::~AliAnalysisTaskSEJPSItoEle()
157 //_________________________________________________________________________________________________
158 void AliAnalysisTaskSEJPSItoEle::Init()
162 if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::Init() \n");
164 gROOT->LoadMacro("ConfigVertexingHF.C");
166 fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()");
171 //_________________________________________________________________________________________________
172 void AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects()
174 // Create the output container
176 if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() \n");
179 Fatal("UserCreateOutputObjects", "This task needs an AOD handler");
183 if(!fNewAOD) fNewAOD = new AliAODEvent();
184 fNewAOD->CreateStdContent();
186 TString filename = "AliAOD.Jpsitoele.root";
187 if (!IsStandardAOD()) filename = "";
189 // Array of HF vertices
190 fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0);
191 fVerticesHFTClArr->SetName("VerticesHF");
192 AddAODBranch("TClonesArray", &fVerticesHFTClArr);
194 // Array of J/psi -> e+e- candidates
195 fJpsiToEleTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0);
196 fJpsiToEleTClArr->SetName("JpsiToEleEvents");
197 AddAODBranch("TClonesArray", &fJpsiToEleTClArr, filename);
199 // Array of like sign candidates
200 fLikeSignTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0);
201 fLikeSignTClArr->SetName("LikeSignEvents");
202 AddAODBranch("TClonesArray", &fLikeSignTClArr, filename);
204 // Array of tracks from J/psi -> e+e- candidates
205 fTracksTClArr = new TClonesArray("AliAODTrack", 0);
206 fTracksTClArr->SetName("Tracks");
207 AddAODBranch("TClonesArray", &fTracksTClArr, filename);
209 fOutput = new TList();
213 fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.);
214 fhDecayTimeMCjpsifromB->Sumw2();
215 fhDecayTimeMCjpsifromB->SetMinimum(0);
216 fOutput->Add(fhDecayTimeMCjpsifromB);
220 fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25);
222 fhInvMass->SetMinimum(0);
223 fOutput->Add(fhInvMass);
225 fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25);
226 fHistMassLS->Sumw2();
227 fHistMassLS->SetMinimum(0);
228 fOutput->Add(fHistMassLS);
230 // TPC signal (only for candidate tracks)
231 fhdEdxTPC = new TH2F("fhdEdxTPC","TPC signal (dE/dx)",100,0.,10.,1000,0.,100.);
232 fhdEdxTPC->SetXTitle("p_{T} [GeV]");
233 fhdEdxTPC->SetYTitle("TPC signal (arb. units)");
234 fOutput->Add(fhdEdxTPC);
236 // Pseudo proper decay time
237 fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.);
238 fhDecayTime->Sumw2();
239 fhDecayTime->SetMinimum(0);
240 fOutput->Add(fhDecayTime);
242 // Pseudo proper decay time
243 fhDecayTimeOut = new TH1F("fhDecayTimeOut", "J/#Psi pseudo proper decay time (output standalone AOD); X [#mu m]; Entries",200,-2000.,2000.);
244 fhDecayTimeOut->Sumw2();
245 fhDecayTimeOut->SetMinimum(0);
246 fOutput->Add(fhDecayTimeOut);
248 // Dictance of closest approach
249 fhDCA = new TH1F("fhDCA", "J/#Psi distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
251 fhDCA->SetMinimum(0);
254 fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.);
256 fHistDCALS->SetMinimum(0);
257 fOutput->Add(fHistDCALS);
260 fhD0 = new TH1F("fhD0", "Impact parameter; d0 [#mu m]; Entries",100,-5000.,5000.);
265 // Product of impact parameters
266 fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.);
268 fhD0D0->SetMinimum(0);
269 fOutput->Add(fhD0D0);
271 fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.);
272 fHistd0d0LS->Sumw2();
273 fHistd0d0LS->SetMinimum(0);
274 fOutput->Add(fHistd0d0LS);
276 // Cosine of the decay angle
277 fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2);
278 fhCosThetaStar->Sumw2();
279 fhCosThetaStar->SetMinimum(0);
280 fOutput->Add(fhCosThetaStar);
282 fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
284 fHistCtsLS->SetMinimum(0);
285 fOutput->Add(fHistCtsLS);
287 fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
288 fHistCtsLSpos->Sumw2();
289 fHistCtsLSpos->SetMinimum(0);
290 fOutput->Add(fHistCtsLSpos);
292 fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.);
293 fHistCtsLSneg->Sumw2();
294 fHistCtsLSneg->SetMinimum(0);
295 fOutput->Add(fHistCtsLSneg);
297 // Cosine of pointing angle
298 fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1);
299 fhCosThetaPointing->Sumw2();
300 fhCosThetaPointing->SetMinimum(0);
301 fOutput->Add(fhCosThetaPointing);
303 fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.);
304 fHistCPtaLS->Sumw2();
305 fHistCPtaLS->SetMinimum(0);
306 fOutput->Add(fHistCPtaLS);
308 // CosThetaStar vs. d0xd0
309 fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000);
310 fhCtsVsD0D0->Sumw2();
311 fhCtsVsD0D0->SetMinimum(0);
312 fOutput->Add(fhCtsVsD0D0);
316 //_________________________________________________________________________________________________
317 void AliAnalysisTaskSEJPSItoEle::UserExec(Option_t */*option*/)
319 // Execute analysis for current event:
321 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
323 TClonesArray *arrayJPSItoEle = 0;
324 TClonesArray *arrayLikeSign = 0;
325 TClonesArray *arrayTracks = 0;
327 if(!aod && AODEvent() && IsStandardAOD()) {
328 // In case there is an AOD handler writing a standard AOD, use the AOD
329 // event in memory rather than the input (ESD) event.
330 aod = dynamic_cast<AliAODEvent*> (AODEvent());
331 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
332 // have to taken from the AOD event hold by the AliAODExtension
333 AliAODHandler* aodHandler = (AliAODHandler*)
334 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
336 if(aodHandler->GetExtensions()) {
337 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
338 AliAODEvent *aodFromExt = ext->GetAOD();
340 // load tracks from AOD event
341 arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks");
342 Int_t totTracks = arrayTracks->GetEntriesFast();
343 if(fDebug>1) printf("Number of tracks === %d \n",totTracks);
345 // load Jpsi candidates from AOD friend
346 arrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle");
347 Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast();
348 if(fDebug>1) printf("Number of J/psi->ee candidates === %d \n",totJPSItoEleCand);
350 // load like sign candidates from AOD friend
351 arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign2Prong");
352 Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast();
353 if(fDebug>1) printf("Number of like sign candidates === %d \n",totLikeSignCand);
357 // load tracks from AOD event
358 arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks");
359 Int_t totTracks = arrayTracks->GetEntriesFast();
360 if(fDebug>1) printf("Number of tracks === %d \n",totTracks);
362 // load Jpsi candidates
363 arrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle");
364 Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast();
365 if(fDebug>1) printf("Number of J/psi->ee candidates === %d \n",totJPSItoEleCand);
367 // load like sign candidates
368 arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong");
369 Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast();
370 if(fDebug>1) printf("Number of like sign candidates === %d \n",totLikeSignCand);
373 fOrigAOD = aod; // copy pointer to the current AliAODEvent in the data member fOrigAOD
377 printf("AliAnalysisTaskSEJPSItoEle::UserExec: Tracks branch not found!\n");
380 if(!arrayJPSItoEle) {
381 printf("AliAnalysisTaskSEJPSItoEle::UserExec: JPSItoEle branch not found!\n");
385 printf("AliAnalysisTaskSEJPSItoEle::UserExec: LikeSign2Prong branch not found!\n");
389 // load MC particles and read MC info (for sim only)
390 TClonesArray* mcArray=0;
392 mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
393 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
394 ReadAODMCInfo(aod,arrayJPSItoEle);
397 // retrieve AOD primary vertex
398 AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
400 Int_t iOutVerticesHF = 0, iOutJPSItoEle = 0, iOutLikeSign = 0, iOutTracks = 0;
402 fVerticesHFTClArr->Delete();
403 iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast();
404 TClonesArray &verticesHFRef = *fVerticesHFTClArr;
406 fJpsiToEleTClArr->Delete();
407 iOutJPSItoEle = fJpsiToEleTClArr->GetEntriesFast();
408 TClonesArray &arrayJpsiToEleRef = *fJpsiToEleTClArr;
410 fLikeSignTClArr->Delete();
411 iOutLikeSign = fLikeSignTClArr->GetEntriesFast();
412 TClonesArray &arrayLikeSignRef = *fLikeSignTClArr;
414 fTracksTClArr->Delete();
415 iOutTracks = fTracksTClArr->GetEntriesFast();
416 //TClonesArray &arrayTrackRef = *fTracksTClArr;
419 // LOOP OVER LIKE SIGN CANDIDATES
422 // Access to the new AOD container of tracks
424 Int_t trkIDtoEntry[100000];
425 AliAODTrack* trackin=0;
426 for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
427 trackin = aod->GetTrack(it);
428 trkIDtoEntry[trackin->GetID()]=it;
431 for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
432 vtx1->AddDaughter(trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack((*(aod->GetTrack(trkIDtoEntry[it])))));
433 trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks);
434 printf("trackin ==== %d \n", trackin);
435 AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin);
437 const AliAODTrack* trackin;
438 Int_t numTracks = arrayTracks->GetEntriesFast();
439 for(Int_t ntracks=0; ntracks<numTracks; ntracks++){
440 trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks);
441 printf("trackin ==== %d \n", trackin);
442 AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin);
446 Int_t nPosPairs=0,nNegPairs=0;
447 Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
448 AliAODRecoDecayHF2Prong* dlsout = 0;
449 //if(fDebug>1) printf("+++\n+++ Number of total like sign pairs (before cuts)---> %d \n+++\n", nLikeSign);
451 for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
452 AliAODRecoDecayHF2Prong *dlsin = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
453 Bool_t unsetvtx=kFALSE;
454 if(!dlsin->GetOwnPrimaryVtx()) {
455 dlsin->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
458 //Int_t okBtoJPSIls=0;
459 //if(dlsin->Pt() < fPtCuts[1] && dlsin->Pt() > fPtCuts[0]){ // apply pt cut only
460 //if(dlsin->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSIls)) {
461 AliAODTrack *trk0 = (AliAODTrack*)dlsin->GetDaughter(0);
462 fHistMassLS->Fill(dlsin->InvMassJPSIee());
463 fHistCPtaLS->Fill(dlsin->CosPointingAngle());
464 fHistd0d0LS->Fill(1e8*dlsin->Prodd0d0());
465 fHistDCALS->Fill(100*dlsin->GetDCA());
467 trk0=aod->GetTrack(trkIDtoEntry[dlsin->GetProngID(0)]);
468 printf("references to standard AOD not available \n");
470 if((trk0->Charge())==1) {
472 fHistCtsLS->Fill(dlsin->CosThetaStar(0,443,11,11));
473 fHistCtsLSpos->Fill(dlsin->CosThetaStar(0,443,11,11));
476 fHistCtsLS->Fill(dlsin->CosThetaStarJPSI());
477 fHistCtsLSneg->Fill(dlsin->CosThetaStarJPSI());
481 // Clone like sign candidate for output AOD
482 dlsout = new(arrayLikeSignRef[iOutLikeSign++]) AliAODRecoDecayHF2Prong(*dlsin);
485 if(unsetvtx) dlsin->UnsetOwnPrimaryVtx();
488 //if(fDebug>1) printf("+++\n+++ N. of positive pairs passing cuts in Event ----- %d \n+++\n", nPosPairs);
489 //if(fDebug>1) printf("+++\n+++ N. of negative pairs passing cuts in Event ----- %d \n+++\n", nNegPairs);
491 fTotPosPairs += nPosPairs;
492 fTotNegPairs += nNegPairs;
495 // LOOP OVER J/psi->e+e- CANDIDATES
498 Int_t nInJPSItoEle = arrayJPSItoEle->GetEntriesFast();
499 //if(fDebug>1) printf("+++\n+++ Number of total like JPSI -> ee candidates (before cuts)---> %d \n+++\n", nInJPSItoEle);
501 //totJPSIin += nInJPSItoEle;
502 Bool_t isOkEle[2] = {kFALSE,kFALSE};
504 for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) {
506 AliAODRecoDecayHF2Prong *dIn = (AliAODRecoDecayHF2Prong*)arrayJPSItoEle->UncheckedAt(iJPSItoEle);
508 if(fOkAODMC) mcLabel = dIn->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
511 Double_t vtxSec[3] = {0.,0.,0.};
512 Double_t vtxPrim[3] = {0.,0.,0.};
513 Double_t vtxSecOut[3] = {0.,0.,0.};
514 Double_t vtxPrimOut[3] = {0.,0.,0.};
515 dIn->GetSecondaryVtx(vtxSec);
516 Bool_t unsetvtx=kFALSE;
517 if(!dIn->GetOwnPrimaryVtx()) {
518 dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
523 //if(dIn->Pt() < fPtCuts[1] && dIn->Pt() > fPtCuts[0]){ // apply pt cut only
524 //if(dIn->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSI)) {
525 if ( fOkAODMC && mcLabel == -1){AliDebug(2,"No MC particle found");} else {
527 //fhInvMass->Fill(dIn->InvMassJPSIee());
528 fhD0->Fill(10000*dIn->ImpParXY());
529 fhD0D0->Fill(1e8*dIn->Prodd0d0());
530 fhCosThetaStar->Fill(dIn->CosThetaStarJPSI());
531 fhCtsVsD0D0->Fill(dIn->CosThetaStarJPSI(),1e8*dIn->Prodd0d0());
532 fhCosThetaPointing->Fill(dIn->CosPointingAngle());
533 fhDCA->Fill(100*dIn->GetDCA());
535 // compute X variable
536 AliAODVertex* primVertex = dIn->GetOwnPrimaryVtx();
537 vtxPrim[0] = primVertex->GetX();
538 vtxPrim[1] = primVertex->GetY();
539 vtxPrim[2] = primVertex->GetZ();
540 Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dIn->Px()) + (vtxSec[1]-vtxPrim[1])*(dIn->Py()))/dIn->Pt();
541 Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dIn->Pt();
542 fhDecayTime->Fill(10000*pseudoProperDecayTime);
544 // retrieve daughter tracks
545 AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0);
546 TBits clusterMap0 = trk0->GetTPCClusterMap();
547 Int_t npoints0 = clusterMap0.CountBits(0);
548 AliAODPid *pid0 = trk0->GetDetPid();
549 if(GetNumberOfSigmas(trk0->P(),pid0->GetTPCsignal(),npoints0,AliPID::kElectron) < 3.) isOkEle[0]=kTRUE;
551 AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1);
552 TBits clusterMap1 = trk1->GetTPCClusterMap();
553 Int_t npoints1 = clusterMap1.CountBits(0);
554 AliAODPid *pid1 = trk1->GetDetPid();
555 if(GetNumberOfSigmas(trk1->P(),pid1->GetTPCsignal(),npoints1,AliPID::kElectron) < 3.) isOkEle[1]=kTRUE;
557 if(isOkEle[0] && isOkEle[1]) {
558 fhInvMass->Fill(dIn->InvMassJPSIee());
559 fhdEdxTPC->Fill(pid0->GetTPCmomentum(),pid0->GetTPCsignal());
560 fhdEdxTPC->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal());
563 //printf("TPC momentum %f\n", pid0->GetTPCmomentum());
564 //printf("TPC signal %f\n", pid0->GetTPCsignal());
566 // clone candidate for output AOD
567 AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++])
568 AliAODVertex(*(dIn->GetSecondaryVtx()));
569 AliAODRecoDecayHF2Prong* dOut = new(arrayJpsiToEleRef[iOutJPSItoEle++])
570 AliAODRecoDecayHF2Prong(*dIn); // copy constructor used
571 dOut->SetSecondaryVtx(v);
572 dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone()));
575 // compute X variable from the cloned candidate in the stand-alone AOD file
576 dOut->GetSecondaryVtx(vtxSecOut);
577 AliAODVertex* primVertexOut = dOut->GetOwnPrimaryVtx();
578 vtxPrimOut[0] = primVertexOut->GetX();
579 vtxPrimOut[1] = primVertexOut->GetY();
580 vtxPrimOut[2] = primVertexOut->GetZ();
581 Double_t lxyOut = ((vtxSecOut[0]-vtxPrimOut[0])*(dOut->Px()) + (vtxSecOut[1]-vtxPrimOut[1])*(dOut->Py()))/dOut->Pt();
582 Double_t pseudoProperDecayTimeOut = lxyOut*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dOut->Pt();
583 fhDecayTimeOut->Fill(10000*pseudoProperDecayTimeOut);
585 //AliAODTrack* trk0 = (AliAODTrack*)dOut->GetDaughter(0);
586 //printf("+++ Pt first daughter = %d \n +++ \n", trk0->Pt());
590 //} // end of JPSItoEle candidates selection according to cuts
592 if(unsetvtx) dIn->UnsetOwnPrimaryVtx();
596 }// end loop on JPSI to ele candidates
598 //printf("+++\n+++ Number of selected J/psi->e+e-: %d\n+++\n",iOutJPSItoEle);
600 //totJPSIout += iOutJPSItoEle;
605 //_________________________________________________________________________________________________
606 void AliAnalysisTaskSEJPSItoEle::Terminate(Option_t */*option*/)
609 // Terminate analysis
612 if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle: Terminate() \n");
614 fOutput = dynamic_cast<TList*> (GetOutputData(1));
616 printf("ERROR: fOutput not available\n");
620 fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs);
622 if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeMCjpsifromB"));
623 fhDecayTime = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTime"));
624 fhDecayTimeOut = dynamic_cast<TH1F*>(fOutput->FindObject("fhDecayTimeOut"));
625 fhInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMass"));
626 fhdEdxTPC = dynamic_cast<TH2F*>(fOutput->FindObject("fhdEdxTPC"));
627 fhD0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0"));
628 fhD0D0 = dynamic_cast<TH1F*>(fOutput->FindObject("fhD0D0"));
629 fhCosThetaStar = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaStar"));
630 fhCosThetaPointing = dynamic_cast<TH1F*>(fOutput->FindObject("fhCosThetaPointing"));
631 fhCtsVsD0D0 = dynamic_cast<TH2F*>(fOutput->FindObject("fhCtsVsD0D0"));
632 fHistMassLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistMassLS"));
633 fHistCtsLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLS"));
634 fHistCtsLSpos = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSpos"));
635 fHistCtsLSneg = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCtsLSneg"));
636 fHistCPtaLS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistCPtaLS"));
637 fHistd0d0LS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistd0d0LS"));
638 fhDCA = dynamic_cast<TH1F*>(fOutput->FindObject("fhDCA"));
639 fHistDCALS = dynamic_cast<TH1F*>(fOutput->FindObject("fHistDCALS"));
641 if(fLsNormalization>0.) {
642 fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries());
643 fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries());
644 fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries());
645 fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries());
646 fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries());
647 fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries());
648 fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries());
651 //printf("Tot JPSI in %d\n", totJPSIin);
652 //printf("Tot JPSI out %d\n", totJPSIout);
656 //_________________________________________________________________________________________________
657 void AliAnalysisTaskSEJPSItoEle::SetPtCuts(const Double_t ptCuts[2])
660 // apply Pt cuts (lower cut, upper cut)
662 for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i];
664 //_________________________________________________________________________________________________
665 void AliAnalysisTaskSEJPSItoEle::SetCutsJPSI(const Double_t cuts[9])
670 for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i];
672 //_________________________________________________________________________________________________
673 void AliAnalysisTaskSEJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray)
676 // Reads MC information if needed (for example in case of parametrized PID)
680 TClonesArray* mcArray =
681 dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
682 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
683 AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast()));
686 AliAODMCHeader* mcHeader =
687 (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
689 printf("AliAnalysisTaskSEJPSItoEle::UserExec: MC AOD header branch not found!\n");
692 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
694 Double_t mcPrimVtx[3];
696 // get MC primary Vtx
697 mcHeader->GetVertex(mcPrimVtx);
699 Int_t nInCandidates = inputArray->GetEntriesFast();
700 printf("Number of Candidates for MC analysis: %d\n",nInCandidates);
702 Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother;
703 Int_t labJPSIdaugh0=0;
704 Int_t labJPSIdaugh1=0;
706 for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) {
708 AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate);
709 Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI
711 Double_t vtxPrim[3] = {0.,0.,0.};
712 Double_t vtxSec[3] = {0.,0.,0.};
713 dd->GetSecondaryVtx(vtxSec);
714 Bool_t unsetvtx=kFALSE;
715 if(!dd->GetOwnPrimaryVtx()) {
716 dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
720 if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only
722 AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0);
723 AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1);
724 lab0 = trk0->GetLabel();
725 lab1 = trk1->GetLabel();
727 AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0);
729 printf("no MC particle\n");
733 while(part0->GetMother()>=0) {
734 labMother=part0->GetMother();
735 part0 = (AliAODMCParticle*)mcArray->At(labMother);
737 printf("no MC mother particle\n");
740 pdgMother = TMath::Abs(part0->GetPdgCode());
741 if(pdgMother==443) {//this for JPSI
742 labJPSIdaugh0=labMother;
747 AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1);
749 printf("no MC particle\n");
753 while(part1->GetMother()>=0) {
754 labMother=part1->GetMother();
755 part1 = (AliAODMCParticle*)mcArray->At(labMother);
757 printf("no MC mother particle\n");
760 pdgMother = TMath::Abs(part1->GetPdgCode());
761 if(pdgMother==443) {//this for JPSI
762 labJPSIdaugh1=labMother;
767 if (mcLabelSec == -1)
769 AliDebug(2,"No MC particle found");
773 if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal
774 AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0);
775 pdgJPSI = partJPSI->GetPdgCode();
776 Int_t pdaugh0 = partJPSI->GetDaughter(0);
777 Int_t pdaugh1 = partJPSI->GetDaughter(1);
778 AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0);
779 AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1);
780 pdg0 = partDaugh0->GetPdgCode();
781 pdg1 = partDaugh1->GetPdgCode();
782 if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee
783 labJPSIMother = partJPSI->GetMother();
784 AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother);
785 pdgJPSIMother = partJPSIMother->GetPdgCode();
787 // chech if JPSI comes from B
788 if( pdgJPSIMother==511 || pdgJPSIMother==521 ||
789 pdgJPSIMother==10511 || pdgJPSIMother==10521 ||
790 pdgJPSIMother==513 || pdgJPSIMother==523 ||
791 pdgJPSIMother==10513 || pdgJPSIMother==10523 ||
792 pdgJPSIMother==20513 || pdgJPSIMother==20523 ||
793 pdgJPSIMother==515 || pdgJPSIMother==525 ||
794 pdgJPSIMother==531 || pdgJPSIMother==10531 ||
795 pdgJPSIMother==533 || pdgJPSIMother==10533 ||
796 pdgJPSIMother==20533 || pdgJPSIMother==535 ||
797 pdgJPSIMother==541 || pdgJPSIMother==10541 ||
798 pdgJPSIMother==543 || pdgJPSIMother==10543 ||
799 pdgJPSIMother==20543 || pdgJPSIMother==545)
800 { //this is for MC JPSI -> ee from B-hadrons
802 fhInvMass->Fill(dd->InvMassJPSIee());
803 fhD0->Fill(10000*dd->ImpParXY());
804 fhD0D0->Fill(1e8*dd->Prodd0d0());
805 fhCosThetaStar->Fill(dd->CosThetaStarJPSI());
806 fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0());
807 fhCosThetaPointing->Fill(dd->CosPointingAngle());
809 // compute X variable
810 AliAODVertex* primVertex = dd->GetOwnPrimaryVtx();
811 vtxPrim[0] = primVertex->GetX();
812 vtxPrim[1] = primVertex->GetY();
813 vtxPrim[2] = primVertex->GetZ();
814 Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt();
815 Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt();
817 fhDecayTime->Fill(10000*pseudoProperDecayTime);
818 // Post the data already here
821 Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt();
822 Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt();
823 fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1);
825 // Post the data already here
828 } //this is for MC JPSI -> ee from B-hadrons
829 } //this is for MC JPSI -> ee
831 } // end of check if JPSI is signal