/************************************************************************** * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ //************************************************************************* // Class AliAnalysisTaskSEJPSItoEle // AliAnalysisTaskSE class to read both J/psi -> e+e- candidates and // like-sign pair candidates and to store them into a specific // stand-alone AOD file for J/psi into dieletrons analysis. // The current task: // // 1) produces several histograms (including invariant mass distributions) // for both unlike sign (US) and like sign (LS) pairs. // // 2) selects only J/Psi to dielectrons candidates and copies it to a // specific AOD file, namely AliAODjpsi.root. The references // to AliAODTrack objects are copied as well. // The specific AliAODjpsi.root is thought as the input file // to the AliAnalysisBtoJPSItoEle.h class, which performs (locally) // the analysis of J/Psi from beauty hadrons decays. // // Author: C.Di Giglio, carmelo.digiglio@ba.infn.it //************************************************************************* class AliAnalysisTaskSE; class AliESDEvent; class AliVEvent; class iostream; class TSystem; class TROOT; #include "TFile.h" #include "TClonesArray.h" #include "TDatabasePDG.h" #include "TROOT.h" #include "TCanvas.h" #include "TBits.h" //********************** #include "AliAODEvent.h" #include "AliAODMCParticle.h" #include "AliAODMCHeader.h" #include "AliAODRecoDecayHF2Prong.h" #include "AliAODHandler.h" #include "AliAnalysisManager.h" #include "AliAnalysisTaskSEJPSItoEle.h" #include "AliAODTrack.h" #include "AliExternalTrackParam.h" ClassImp(AliAnalysisTaskSEJPSItoEle) //______________________________________________________________________________ AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle() : AliAnalysisTaskSE(), fOutput(0), fhDecayTimeMCjpsifromB(0), fhDecayTime(0), fhDecayTimeOut(0), fhInvMass(0), fhdEdxTPC(0), fhD0(0), fhD0D0(0), fhCosThetaStar(0), fhCosThetaPointing(0), fhDCA(0), fhCtsVsD0D0(0), fHistMassLS(0), fHistCtsLS(0), fHistCtsLSpos(0), fHistCtsLSneg(0), fHistCPtaLS(0), fHistd0d0LS(0), fHistDCALS(0), fVHF(0), fTotPosPairs(0), fTotNegPairs(0), fLsNormalization(1.), fOkAODMC(kFALSE), fOkLikeSign(kFALSE), fVerticesHFTClArr(0), fJpsiToEleTClArr(0), fLikeSignTClArr(0), fTracksTClArr(0), fOrigAOD(0), fNewAOD(0) { // default constructor } //_________________________________________________________________________________________________ AliAnalysisTaskSEJPSItoEle::AliAnalysisTaskSEJPSItoEle(const char* name) : AliAnalysisTaskSE(name), fOutput(0), fhDecayTimeMCjpsifromB(0), fhDecayTime(0), fhDecayTimeOut(0), fhInvMass(0), fhdEdxTPC(0), fhD0(0), fhD0D0(0), fhCosThetaStar(0), fhCosThetaPointing(0), fhDCA(0), fhCtsVsD0D0(0), fHistMassLS(0), fHistCtsLS(0), fHistCtsLSpos(0), fHistCtsLSneg(0), fHistCPtaLS(0), fHistd0d0LS(0), fHistDCALS(0), fVHF(0), fTotPosPairs(0), fTotNegPairs(0), fLsNormalization(1.), fOkAODMC(kFALSE), fOkLikeSign(kFALSE), fVerticesHFTClArr(0), fJpsiToEleTClArr(0), fLikeSignTClArr(0), fTracksTClArr(0), fOrigAOD(0), fNewAOD(0) { // standard constructor // Input slot #0 works with an Ntuple DefineInput(0, TChain::Class()); // Output slot #0 writes into a TTree container // Output slot #1 writes into a TList container DefineOutput(0, TTree::Class()); DefineOutput(1,TList::Class()); //My private output } //_________________________________________________________________________________________________ AliAnalysisTaskSEJPSItoEle::~AliAnalysisTaskSEJPSItoEle() { // destructor if (fOutput) { delete fOutput; fOutput = 0; } if (fVHF) { delete fVHF; fVHF = 0; } } //_________________________________________________________________________________________________ void AliAnalysisTaskSEJPSItoEle::Init() { // Initialization if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::Init() \n"); gROOT->LoadMacro("ConfigVertexingHF.C"); fVHF = (AliAnalysisVertexingHF*)gROOT->ProcessLine("ConfigVertexingHF()"); fVHF->PrintStatus(); return; } //_________________________________________________________________________________________________ void AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() { // Create the output container if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle::UserCreateOutputObjects() \n"); if (!AODEvent()) { Fatal("UserCreateOutputObjects", "This task needs an AOD handler"); return; } if(!fNewAOD) fNewAOD = new AliAODEvent(); fNewAOD->CreateStdContent(); TString filename = "AliAOD.Jpsitoele.root"; if (!IsStandardAOD()) filename = ""; // Array of HF vertices fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0); fVerticesHFTClArr->SetName("VerticesHF"); AddAODBranch("TClonesArray", &fVerticesHFTClArr); // Array of J/psi -> e+e- candidates fJpsiToEleTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0); fJpsiToEleTClArr->SetName("JpsiToEleEvents"); AddAODBranch("TClonesArray", &fJpsiToEleTClArr, filename); // Array of like sign candidates fLikeSignTClArr = new TClonesArray("AliAODRecoDecayHF2Prong",0); fLikeSignTClArr->SetName("LikeSignEvents"); AddAODBranch("TClonesArray", &fLikeSignTClArr, filename); // Array of tracks from J/psi -> e+e- candidates fTracksTClArr = new TClonesArray("AliAODTrack", 0); fTracksTClArr->SetName("Tracks"); AddAODBranch("TClonesArray", &fTracksTClArr, filename); fOutput = new TList(); fOutput->SetOwner(); if(fOkAODMC){ fhDecayTimeMCjpsifromB = new TH1F("fhDecayTimeMCjpsifromB", "Secondary J/#Psi Monte Carlo pseudo proper decay time; X [#mu m]; Entries",100,-4000.,4000.); fhDecayTimeMCjpsifromB->Sumw2(); fhDecayTimeMCjpsifromB->SetMinimum(0); fOutput->Add(fhDecayTimeMCjpsifromB); } // Invariant mass fhInvMass = new TH1F("fhInvMass", "J/#Psi invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25); fhInvMass->Sumw2(); fhInvMass->SetMinimum(0); fOutput->Add(fhInvMass); fHistMassLS = new TH1F("fHistMassLS", "Like sign pairs invariant mass; Inv. mass [GeV]; Entries/5 MeV",100,2.75,3.25); fHistMassLS->Sumw2(); fHistMassLS->SetMinimum(0); fOutput->Add(fHistMassLS); // TPC signal (only for candidate tracks) fhdEdxTPC = new TH2F("fhdEdxTPC","TPC signal (dE/dx)",100,0.,10.,1000,0.,100.); fhdEdxTPC->SetXTitle("p_{T} [GeV]"); fhdEdxTPC->SetYTitle("TPC signal (arb. units)"); fOutput->Add(fhdEdxTPC); // Pseudo proper decay time fhDecayTime = new TH1F("fhDecayTime", "J/#Psi pseudo proper decay time; X [#mu m]; Entries",200,-2000.,2000.); fhDecayTime->Sumw2(); fhDecayTime->SetMinimum(0); fOutput->Add(fhDecayTime); // Pseudo proper decay time fhDecayTimeOut = new TH1F("fhDecayTimeOut", "J/#Psi pseudo proper decay time (output standalone AOD); X [#mu m]; Entries",200,-2000.,2000.); fhDecayTimeOut->Sumw2(); fhDecayTimeOut->SetMinimum(0); fOutput->Add(fhDecayTimeOut); // Dictance of closest approach fhDCA = new TH1F("fhDCA", "J/#Psi distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.); fhDCA->Sumw2(); fhDCA->SetMinimum(0); fOutput->Add(fhDCA); fHistDCALS = new TH1F("fHistDCALS", "Like sign pairs distance of closest approach; dca [10^{2}#mu m]; Entries",100,0.,5.); fHistDCALS->Sumw2(); fHistDCALS->SetMinimum(0); fOutput->Add(fHistDCALS); // Impact parameter fhD0 = new TH1F("fhD0", "Impact parameter; d0 [#mu m]; Entries",100,-5000.,5000.); fhD0->Sumw2(); fhD0->SetMinimum(0); fOutput->Add(fhD0); // Product of impact parameters fhD0D0 = new TH1F("fhD0D0", "Product of impact parameters; D0D0 [#mu m^2]; Entries",100,-100000.,100000.); fhD0D0->Sumw2(); fhD0D0->SetMinimum(0); fOutput->Add(fhD0D0); fHistd0d0LS = new TH1F("fHistd0d0LS", "Like sign pairs product of impact parameters; d0xd0 [#mu m^{2}]; Entries",200,-100000.,100000.); fHistd0d0LS->Sumw2(); fHistd0d0LS->SetMinimum(0); fOutput->Add(fHistd0d0LS); // Cosine of the decay angle fhCosThetaStar = new TH1F("fhCosThetaStar", "Cosine of decay angle; Cosine Theta Star; Entries",50,-1.2,1.2); fhCosThetaStar->Sumw2(); fhCosThetaStar->SetMinimum(0); fOutput->Add(fhCosThetaStar); fHistCtsLS = new TH1F("fHistCtsLS", "Like sign pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); fHistCtsLS->Sumw2(); fHistCtsLS->SetMinimum(0); fOutput->Add(fHistCtsLS); fHistCtsLSpos = new TH1F("fHistCtsLSpos", "Like sign ++ pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); fHistCtsLSpos->Sumw2(); fHistCtsLSpos->SetMinimum(0); fOutput->Add(fHistCtsLSpos); fHistCtsLSneg = new TH1F("fHistCtsLSneg", "Like sign -- pairs cosine of decay angle; Cos#Theta*; Entries",200,-1.,1.); fHistCtsLSneg->Sumw2(); fHistCtsLSneg->SetMinimum(0); fOutput->Add(fHistCtsLSneg); // Cosine of pointing angle fhCosThetaPointing = new TH1F("fhCosThetaPointing", "Cosine of pointing angle; Cosine Pointing Angle; Entries",100,-1,1); fhCosThetaPointing->Sumw2(); fhCosThetaPointing->SetMinimum(0); fOutput->Add(fhCosThetaPointing); fHistCPtaLS = new TH1F("fHistCPtaLS", "Like sign pairs cosine of pointing angle; Cos#Theta_{point}; Entries",200,-1.,1.); fHistCPtaLS->Sumw2(); fHistCPtaLS->SetMinimum(0); fOutput->Add(fHistCPtaLS); // CosThetaStar vs. d0xd0 fhCtsVsD0D0 = new TH2F("fhCtsVsD0D0", "Cosine theta star Vs. D0; Cosine theta star; D0D0",50,-1,1,100,-100000,100000); fhCtsVsD0D0->Sumw2(); fhCtsVsD0D0->SetMinimum(0); fOutput->Add(fhCtsVsD0D0); return; } //_________________________________________________________________________________________________ void AliAnalysisTaskSEJPSItoEle::UserExec(Option_t */*option*/) { // Execute analysis for current event: AliAODEvent *aod = dynamic_cast (InputEvent()); TClonesArray *arrayJPSItoEle = 0; TClonesArray *arrayLikeSign = 0; TClonesArray *arrayTracks = 0; if(!aod && AODEvent() && IsStandardAOD()) { // In case there is an AOD handler writing a standard AOD, use the AOD // event in memory rather than the input (ESD) event. aod = dynamic_cast (AODEvent()); // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root) // have to taken from the AOD event hold by the AliAODExtension AliAODHandler* aodHandler = (AliAODHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()); if(aodHandler->GetExtensions()) { AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root"); AliAODEvent *aodFromExt = ext->GetAOD(); // load tracks from AOD event arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks"); Int_t totTracks = arrayTracks->GetEntriesFast(); if(fDebug>1) printf("Number of tracks === %d \n",totTracks); // load Jpsi candidates from AOD friend arrayJPSItoEle=(TClonesArray*)aodFromExt->GetList()->FindObject("JPSItoEle"); Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast(); if(fDebug>1) printf("Number of J/psi->ee candidates === %d \n",totJPSItoEleCand); // load like sign candidates from AOD friend arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign2Prong"); Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast(); if(fDebug>1) printf("Number of like sign candidates === %d \n",totLikeSignCand); } } else { // load tracks from AOD event arrayTracks=(TClonesArray*)aod->GetList()->FindObject("tracks"); Int_t totTracks = arrayTracks->GetEntriesFast(); if(fDebug>1) printf("Number of tracks === %d \n",totTracks); // load Jpsi candidates arrayJPSItoEle=(TClonesArray*)aod->GetList()->FindObject("JPSItoEle"); Int_t totJPSItoEleCand = arrayJPSItoEle->GetEntriesFast(); if(fDebug>1) printf("Number of J/psi->ee candidates === %d \n",totJPSItoEleCand); // load like sign candidates arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign2Prong"); Int_t totLikeSignCand = arrayLikeSign->GetEntriesFast(); if(fDebug>1) printf("Number of like sign candidates === %d \n",totLikeSignCand); } fOrigAOD = aod; // copy pointer to the current AliAODEvent in the data member fOrigAOD if (!aod) return; if(!arrayTracks) { printf("AliAnalysisTaskSEJPSItoEle::UserExec: Tracks branch not found!\n"); return; } if(!arrayJPSItoEle) { printf("AliAnalysisTaskSEJPSItoEle::UserExec: JPSItoEle branch not found!\n"); return; } if(!arrayLikeSign) { printf("AliAnalysisTaskSEJPSItoEle::UserExec: LikeSign2Prong branch not found!\n"); return; } // fix for temporary bug in ESDfilter // the AODs with null vertex pointer didn't pass the PhysSel if(!aod->GetPrimaryVertex()) return; // load MC particles and read MC info (for sim only) TClonesArray* mcArray=0; if(fOkAODMC){ mcArray = dynamic_cast(aod->FindListObject(AliAODMCParticle::StdBranchName())); if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); ReadAODMCInfo(aod,arrayJPSItoEle); } // retrieve AOD primary vertex AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex(); Int_t iOutVerticesHF = 0, iOutJPSItoEle = 0, iOutLikeSign = 0, iOutTracks = 0; fVerticesHFTClArr->Delete(); iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast(); TClonesArray &verticesHFRef = *fVerticesHFTClArr; fJpsiToEleTClArr->Delete(); iOutJPSItoEle = fJpsiToEleTClArr->GetEntriesFast(); TClonesArray &arrayJpsiToEleRef = *fJpsiToEleTClArr; fLikeSignTClArr->Delete(); iOutLikeSign = fLikeSignTClArr->GetEntriesFast(); TClonesArray &arrayLikeSignRef = *fLikeSignTClArr; fTracksTClArr->Delete(); iOutTracks = fTracksTClArr->GetEntriesFast(); //TClonesArray &arrayTrackRef = *fTracksTClArr; // // LOOP OVER LIKE SIGN CANDIDATES // // Access to the new AOD container of tracks Int_t trkIDtoEntry[100000]; AliAODTrack* trackin=0; for(Int_t it=0;itGetNumberOfTracks();it++) { trackin = aod->GetTrack(it); trkIDtoEntry[trackin->GetID()]=it; } /* for(Int_t it=0;itGetNumberOfTracks();it++) { vtx1->AddDaughter(trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack((*(aod->GetTrack(trkIDtoEntry[it]))))); trackin = (AliAODTrack*)arrayTracks->UncheckedAt(ntracks); printf("trackin ==== %d \n", trackin); AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin); } const AliAODTrack* trackin; Int_t numTracks = arrayTracks->GetEntriesFast(); for(Int_t ntracks=0; ntracksUncheckedAt(ntracks); printf("trackin ==== %d \n", trackin); AliAODTrack* trackout = new(arrayTrackRef[iOutTracks++]) AliAODTrack(*trackin); } */ Int_t nPosPairs=0,nNegPairs=0; Int_t nLikeSign = arrayLikeSign->GetEntriesFast(); AliAODRecoDecayHF2Prong* dlsout = 0; //if(fDebug>1) printf("+++\n+++ Number of total like sign pairs (before cuts)---> %d \n+++\n", nLikeSign); for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) { AliAODRecoDecayHF2Prong *dlsin = (AliAODRecoDecayHF2Prong*)arrayLikeSign->UncheckedAt(iLikeSign); Bool_t unsetvtx=kFALSE; if(!dlsin->GetOwnPrimaryVtx()) { dlsin->SetOwnPrimaryVtx(vtx1); // needed to compute all variables unsetvtx=kTRUE; } //Int_t okBtoJPSIls=0; //if(dlsin->Pt() < fPtCuts[1] && dlsin->Pt() > fPtCuts[0]){ // apply pt cut only //if(dlsin->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSIls)) { AliAODTrack *trk0 = (AliAODTrack*)dlsin->GetDaughter(0); fHistMassLS->Fill(dlsin->InvMassJPSIee()); fHistCPtaLS->Fill(dlsin->CosPointingAngle()); fHistd0d0LS->Fill(1e8*dlsin->Prodd0d0()); fHistDCALS->Fill(100*dlsin->GetDCA()); if(!trk0) { trk0=aod->GetTrack(trkIDtoEntry[dlsin->GetProngID(0)]); printf("references to standard AOD not available \n"); } if((trk0->Charge())==1) { nPosPairs++; fHistCtsLS->Fill(dlsin->CosThetaStar(0,443,11,11)); fHistCtsLSpos->Fill(dlsin->CosThetaStar(0,443,11,11)); } else { nNegPairs++; fHistCtsLS->Fill(dlsin->CosThetaStarJPSI()); fHistCtsLSneg->Fill(dlsin->CosThetaStarJPSI()); } PostData(1,fOutput); // Clone like sign candidate for output AOD dlsout = new(arrayLikeSignRef[iOutLikeSign++]) AliAODRecoDecayHF2Prong(*dlsin); //} if(unsetvtx) dlsin->UnsetOwnPrimaryVtx(); } //if(fDebug>1) printf("+++\n+++ N. of positive pairs passing cuts in Event ----- %d \n+++\n", nPosPairs); //if(fDebug>1) printf("+++\n+++ N. of negative pairs passing cuts in Event ----- %d \n+++\n", nNegPairs); fTotPosPairs += nPosPairs; fTotNegPairs += nNegPairs; // // LOOP OVER J/psi->e+e- CANDIDATES // Int_t nInJPSItoEle = arrayJPSItoEle->GetEntriesFast(); //if(fDebug>1) printf("+++\n+++ Number of total like JPSI -> ee candidates (before cuts)---> %d \n+++\n", nInJPSItoEle); //totJPSIin += nInJPSItoEle; Bool_t isOkEle[2] = {kFALSE,kFALSE}; for (Int_t iJPSItoEle = 0; iJPSItoEle < nInJPSItoEle; iJPSItoEle++) { AliAODRecoDecayHF2Prong *dIn = (AliAODRecoDecayHF2Prong*)arrayJPSItoEle->UncheckedAt(iJPSItoEle); Int_t mcLabel = 0; if(fOkAODMC) mcLabel = dIn->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI //Secondary vertex Double_t vtxSec[3] = {0.,0.,0.}; Double_t vtxPrim[3] = {0.,0.,0.}; Double_t vtxSecOut[3] = {0.,0.,0.}; Double_t vtxPrimOut[3] = {0.,0.,0.}; dIn->GetSecondaryVtx(vtxSec); Bool_t unsetvtx=kFALSE; if(!dIn->GetOwnPrimaryVtx()) { dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables unsetvtx=kTRUE; } //Int_t okBtoJPSI=0; //if(dIn->Pt() < fPtCuts[1] && dIn->Pt() > fPtCuts[0]){ // apply pt cut only //if(dIn->SelectBtoJPSI(fVHF->GetBtoJPSICuts(),okBtoJPSI)) { if ( fOkAODMC && mcLabel == -1){AliDebug(2,"No MC particle found");} else { //fhInvMass->Fill(dIn->InvMassJPSIee()); fhD0->Fill(10000*dIn->ImpParXY()); fhD0D0->Fill(1e8*dIn->Prodd0d0()); fhCosThetaStar->Fill(dIn->CosThetaStarJPSI()); fhCtsVsD0D0->Fill(dIn->CosThetaStarJPSI(),1e8*dIn->Prodd0d0()); fhCosThetaPointing->Fill(dIn->CosPointingAngle()); fhDCA->Fill(100*dIn->GetDCA()); // compute X variable AliAODVertex* primVertex = dIn->GetOwnPrimaryVtx(); vtxPrim[0] = primVertex->GetX(); vtxPrim[1] = primVertex->GetY(); vtxPrim[2] = primVertex->GetZ(); Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dIn->Px()) + (vtxSec[1]-vtxPrim[1])*(dIn->Py()))/dIn->Pt(); Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dIn->Pt(); fhDecayTime->Fill(10000*pseudoProperDecayTime); // retrieve daughter tracks AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0); TBits clusterMap0 = trk0->GetTPCClusterMap(); Int_t npoints0 = clusterMap0.CountBits(0); AliAODPid *pid0 = trk0->GetDetPid(); if(GetNumberOfSigmas(trk0->P(),pid0->GetTPCsignal(),npoints0,AliPID::kElectron) < 3.) isOkEle[0]=kTRUE; AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1); TBits clusterMap1 = trk1->GetTPCClusterMap(); Int_t npoints1 = clusterMap1.CountBits(0); AliAODPid *pid1 = trk1->GetDetPid(); if(GetNumberOfSigmas(trk1->P(),pid1->GetTPCsignal(),npoints1,AliPID::kElectron) < 3.) isOkEle[1]=kTRUE; if(isOkEle[0] && isOkEle[1]) { fhInvMass->Fill(dIn->InvMassJPSIee()); fhdEdxTPC->Fill(pid0->GetTPCmomentum(),pid0->GetTPCsignal()); fhdEdxTPC->Fill(pid1->GetTPCmomentum(),pid1->GetTPCsignal()); } //printf("TPC momentum %f\n", pid0->GetTPCmomentum()); //printf("TPC signal %f\n", pid0->GetTPCsignal()); // clone candidate for output AOD AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++]) AliAODVertex(*(dIn->GetSecondaryVtx())); AliAODRecoDecayHF2Prong* dOut = new(arrayJpsiToEleRef[iOutJPSItoEle++]) AliAODRecoDecayHF2Prong(*dIn); // copy constructor used dOut->SetSecondaryVtx(v); dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone())); v->SetParent(dOut); // compute X variable from the cloned candidate in the stand-alone AOD file dOut->GetSecondaryVtx(vtxSecOut); AliAODVertex* primVertexOut = dOut->GetOwnPrimaryVtx(); vtxPrimOut[0] = primVertexOut->GetX(); vtxPrimOut[1] = primVertexOut->GetY(); vtxPrimOut[2] = primVertexOut->GetZ(); Double_t lxyOut = ((vtxSecOut[0]-vtxPrimOut[0])*(dOut->Px()) + (vtxSecOut[1]-vtxPrimOut[1])*(dOut->Py()))/dOut->Pt(); Double_t pseudoProperDecayTimeOut = lxyOut*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dOut->Pt(); fhDecayTimeOut->Fill(10000*pseudoProperDecayTimeOut); //AliAODTrack* trk0 = (AliAODTrack*)dOut->GetDaughter(0); //printf("+++ Pt first daughter = %d \n +++ \n", trk0->Pt()); } //} // end of JPSItoEle candidates selection according to cuts if(unsetvtx) dIn->UnsetOwnPrimaryVtx(); PostData(0,fOutput); }// end loop on JPSI to ele candidates //printf("+++\n+++ Number of selected J/psi->e+e-: %d\n+++\n",iOutJPSItoEle); //totJPSIout += iOutJPSItoEle; return; } //_________________________________________________________________________________________________ void AliAnalysisTaskSEJPSItoEle::Terminate(Option_t */*option*/) { // // Terminate analysis // if(fDebug > 1) printf("AliAnalysisTaskSEJPSItoEle: Terminate() \n"); fOutput = dynamic_cast (GetOutputData(1)); if (!fOutput) { printf("ERROR: fOutput not available\n"); return; } fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs*fTotNegPairs); if(fOkAODMC) fhDecayTimeMCjpsifromB = dynamic_cast(fOutput->FindObject("fhDecayTimeMCjpsifromB")); fhDecayTime = dynamic_cast(fOutput->FindObject("fhDecayTime")); fhDecayTimeOut = dynamic_cast(fOutput->FindObject("fhDecayTimeOut")); fhInvMass = dynamic_cast(fOutput->FindObject("fhInvMass")); fhdEdxTPC = dynamic_cast(fOutput->FindObject("fhdEdxTPC")); fhD0 = dynamic_cast(fOutput->FindObject("fhD0")); fhD0D0 = dynamic_cast(fOutput->FindObject("fhD0D0")); fhCosThetaStar = dynamic_cast(fOutput->FindObject("fhCosThetaStar")); fhCosThetaPointing = dynamic_cast(fOutput->FindObject("fhCosThetaPointing")); fhCtsVsD0D0 = dynamic_cast(fOutput->FindObject("fhCtsVsD0D0")); fHistMassLS = dynamic_cast(fOutput->FindObject("fHistMassLS")); fHistCtsLS = dynamic_cast(fOutput->FindObject("fHistCtsLS")); fHistCtsLSpos = dynamic_cast(fOutput->FindObject("fHistCtsLSpos")); fHistCtsLSneg = dynamic_cast(fOutput->FindObject("fHistCtsLSneg")); fHistCPtaLS = dynamic_cast(fOutput->FindObject("fHistCPtaLS")); fHistd0d0LS = dynamic_cast(fOutput->FindObject("fHistd0d0LS")); fhDCA = dynamic_cast(fOutput->FindObject("fhDCA")); fHistDCALS = dynamic_cast(fOutput->FindObject("fHistDCALS")); if(fLsNormalization>0.) { fHistMassLS->Scale((1/fLsNormalization)*fHistMassLS->GetEntries()); fHistCtsLS->Scale((1/fLsNormalization)*fHistCtsLS->GetEntries()); fHistCtsLSpos->Scale((1/fLsNormalization)*fHistCtsLSpos->GetEntries()); fHistCtsLSneg->Scale((1/fLsNormalization)*fHistCtsLSneg->GetEntries()); fHistCPtaLS->Scale((1/fLsNormalization)*fHistCPtaLS->GetEntries()); fHistd0d0LS->Scale((1/fLsNormalization)*fHistd0d0LS->GetEntries()); fHistDCALS->Scale((1/fLsNormalization)*fHistDCALS->GetEntries()); } //printf("Tot JPSI in %d\n", totJPSIin); //printf("Tot JPSI out %d\n", totJPSIout); return; } //_________________________________________________________________________________________________ void AliAnalysisTaskSEJPSItoEle::SetPtCuts(const Double_t ptCuts[2]) { // // apply Pt cuts (lower cut, upper cut) // for(Int_t i = 0; i < 2; i++) fPtCuts[i] = ptCuts[i]; } //_________________________________________________________________________________________________ void AliAnalysisTaskSEJPSItoEle::SetCutsJPSI(const Double_t cuts[9]) { // // apply JPSI cuts // for(Int_t i = 0; i < 9; i++) fCuts[i] = cuts[i]; } //_________________________________________________________________________________________________ void AliAnalysisTaskSEJPSItoEle::ReadAODMCInfo(AliAODEvent* aodEvent, const TClonesArray* inputArray) { // // Reads MC information if needed (for example in case of parametrized PID) // // load MC particles TClonesArray* mcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); AliInfo(Form("MC particles found in mcArray ---> %d ",mcArray->GetEntriesFast())); // load MC header AliAODMCHeader* mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()); if(!mcHeader){ printf("AliAnalysisTaskSEJPSItoEle::UserExec: MC AOD header branch not found!\n"); } AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex(); Double_t mcPrimVtx[3]; // get MC primary Vtx mcHeader->GetVertex(mcPrimVtx); Int_t nInCandidates = inputArray->GetEntriesFast(); printf("Number of Candidates for MC analysis: %d\n",nInCandidates); Int_t lab0, lab1, pdgMother, labMother, pdgJPSI, pdg0, pdg1, labJPSIMother, pdgJPSIMother; Int_t labJPSIdaugh0=0; Int_t labJPSIdaugh1=0; for (Int_t iCandidate = 0; iCandidate < nInCandidates; iCandidate++) { AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iCandidate); Int_t mcLabelSec = dd->MatchToMC(443,mcArray) ; // select only reco JPSI that are true JPSI Double_t vtxPrim[3] = {0.,0.,0.}; Double_t vtxSec[3] = {0.,0.,0.}; dd->GetSecondaryVtx(vtxSec); Bool_t unsetvtx=kFALSE; if(!dd->GetOwnPrimaryVtx()) { dd->SetOwnPrimaryVtx(vtx1); // needed to compute all variables unsetvtx=kTRUE; } if(dd->Pt() < fPtCuts[1] && dd->Pt() > fPtCuts[0]){ // apply pt cut only AliAODTrack *trk0 = (AliAODTrack*)dd->GetDaughter(0); AliAODTrack *trk1 = (AliAODTrack*)dd->GetDaughter(1); lab0 = trk0->GetLabel(); lab1 = trk1->GetLabel(); AliAODMCParticle* part0 = (AliAODMCParticle*)mcArray->At(lab0); if(!part0) { printf("no MC particle\n"); continue; } while(part0->GetMother()>=0) { labMother=part0->GetMother(); part0 = (AliAODMCParticle*)mcArray->At(labMother); if(!part0) { printf("no MC mother particle\n"); break; } pdgMother = TMath::Abs(part0->GetPdgCode()); if(pdgMother==443) {//this for JPSI labJPSIdaugh0=labMother; break; } } AliAODMCParticle* part1 = (AliAODMCParticle*)mcArray->At(lab1); if(!part1) { printf("no MC particle\n"); continue; } while(part1->GetMother()>=0) { labMother=part1->GetMother(); part1 = (AliAODMCParticle*)mcArray->At(labMother); if(!part1) { printf("no MC mother particle\n"); break; } pdgMother = TMath::Abs(part1->GetPdgCode()); if(pdgMother==443) {//this for JPSI labJPSIdaugh1=labMother; break; } } if (mcLabelSec == -1) { AliDebug(2,"No MC particle found"); } else { if(labJPSIdaugh0>=0 && labJPSIdaugh1>=0 && labJPSIdaugh0==labJPSIdaugh1) { // check if JPSI is signal AliAODMCParticle* partJPSI = (AliAODMCParticle*)mcArray->At(labJPSIdaugh0); pdgJPSI = partJPSI->GetPdgCode(); Int_t pdaugh0 = partJPSI->GetDaughter(0); Int_t pdaugh1 = partJPSI->GetDaughter(1); AliAODMCParticle* partDaugh0 = (AliAODMCParticle*)mcArray->At(pdaugh0); AliAODMCParticle* partDaugh1 = (AliAODMCParticle*)mcArray->At(pdaugh1); pdg0 = partDaugh0->GetPdgCode(); pdg1 = partDaugh1->GetPdgCode(); if(TMath::Abs(pdg0) == 11 && TMath::Abs(pdg1) == 11){ // this is for MC JPSI -> ee labJPSIMother = partJPSI->GetMother(); AliAODMCParticle* partJPSIMother = (AliAODMCParticle*)mcArray->At(labJPSIMother); pdgJPSIMother = partJPSIMother->GetPdgCode(); // chech if JPSI comes from B if( pdgJPSIMother==511 || pdgJPSIMother==521 || pdgJPSIMother==10511 || pdgJPSIMother==10521 || pdgJPSIMother==513 || pdgJPSIMother==523 || pdgJPSIMother==10513 || pdgJPSIMother==10523 || pdgJPSIMother==20513 || pdgJPSIMother==20523 || pdgJPSIMother==515 || pdgJPSIMother==525 || pdgJPSIMother==531 || pdgJPSIMother==10531 || pdgJPSIMother==533 || pdgJPSIMother==10533 || pdgJPSIMother==20533 || pdgJPSIMother==535 || pdgJPSIMother==541 || pdgJPSIMother==10541 || pdgJPSIMother==543 || pdgJPSIMother==10543 || pdgJPSIMother==20543 || pdgJPSIMother==545) { //this is for MC JPSI -> ee from B-hadrons fhInvMass->Fill(dd->InvMassJPSIee()); fhD0->Fill(10000*dd->ImpParXY()); fhD0D0->Fill(1e8*dd->Prodd0d0()); fhCosThetaStar->Fill(dd->CosThetaStarJPSI()); fhCtsVsD0D0->Fill(dd->CosThetaStarJPSI(),1e8*dd->Prodd0d0()); fhCosThetaPointing->Fill(dd->CosPointingAngle()); // compute X variable AliAODVertex* primVertex = dd->GetOwnPrimaryVtx(); vtxPrim[0] = primVertex->GetX(); vtxPrim[1] = primVertex->GetY(); vtxPrim[2] = primVertex->GetZ(); Double_t lxy = ((vtxSec[0]-vtxPrim[0])*(dd->Px()) + (vtxSec[1]-vtxPrim[1])*(dd->Py()))/dd->Pt(); Double_t pseudoProperDecayTime = lxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/dd->Pt(); fhDecayTime->Fill(10000*pseudoProperDecayTime); // Post the data already here PostData(1,fOutput); Double_t mcLxy = ((partJPSI->Xv()-mcPrimVtx[0])*(partJPSI->Px()) + (partJPSI->Yv()-mcPrimVtx[1])*(partJPSI->Py()))/partJPSI->Pt(); Double_t mcPseudoProperDecayTime = mcLxy*(TDatabasePDG::Instance()->GetParticle(443)->Mass())/partJPSI->Pt(); fhDecayTimeMCjpsifromB->Fill(10000*mcPseudoProperDecayTime,1); // Post the data already here PostData(1,fOutput); } //this is for MC JPSI -> ee from B-hadrons } //this is for MC JPSI -> ee } } // end of check if JPSI is signal } } }