/************************************************************************** * Copyright(c) 1998-1999, 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 for heavy-flavour electron with EMCal triggered events // Author: Shingo Sakai #include "TChain.h" #include "TTree.h" #include "TH2F.h" #include "TF1.h" #include "TMath.h" #include "TCanvas.h" #include "THnSparse.h" #include "TLorentzVector.h" #include "TString.h" #include "TFile.h" #include "TGraphErrors.h" #include "TDatabasePDG.h" #include "AliAnalysisTask.h" #include "AliAnalysisManager.h" #include "AliESDEvent.h" #include "AliESDHandler.h" #include "AliAODEvent.h" #include "AliAODHandler.h" #include "AliMCEventHandler.h" #include "AliMCEvent.h" #include "AliMCParticle.h" #include "AliAnalysisTaskHFECal.h" #include "TGeoGlobalMagField.h" #include "AliLog.h" #include "AliAnalysisTaskSE.h" #include "TRefArray.h" #include "TVector.h" #include "AliESDInputHandler.h" #include "AliESDpid.h" #include "AliESDtrackCuts.h" #include "AliPhysicsSelection.h" #include "AliESDCaloCluster.h" #include "AliAODCaloCluster.h" #include "AliEMCALRecoUtils.h" #include "AliEMCALGeometry.h" #include "AliGeomManager.h" #include "stdio.h" #include "TGeoManager.h" #include "iostream" #include "fstream" #include "AliEMCALTrack.h" #include "AliMagF.h" #include "AliKFParticle.h" #include "AliKFVertex.h" #include "AliPID.h" #include "AliPIDResponse.h" #include "AliHFEcontainer.h" #include "AliHFEcuts.h" #include "AliHFEpid.h" #include "AliHFEpidBase.h" #include "AliHFEpidQAmanager.h" #include "AliHFEtools.h" #include "AliCFContainer.h" #include "AliCFManager.h" #include "AliStack.h" #include "AliCentrality.h" using namespace std; ClassImp(AliAnalysisTaskHFECal) //________________________________________________________________________ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) : AliAnalysisTaskSE(name) ,fESD(0) ,fMC(0) ,stack(0) ,fGeom(0) ,fOutputList(0) ,fqahist(0) ,fTrackCuts(0) ,fCuts(0) ,fIdentifiedAsOutInz(kFALSE) ,fPassTheEventCut(kFALSE) ,fRejectKinkMother(kFALSE) ,fmcData(kFALSE) ,fVz(0.0) ,fCFM(0) ,fPID(0) ,fPIDqa(0) ,fOpeningAngleCut(0.1) ,fMimpTassCut(0.5) ,fMimNsigassCut(-3) ,fInvmassCut(0) // no mass ,fSetMassConstraint(kTRUE) ,fSetMassWidthCut(kFALSE) ,fSetMassNonlinear(kFALSE) ,fSetKFpart(kTRUE) ,fNoEvents(0) ,fEMCAccE(0) ,hEMCAccE(0) ,fTrkpt(0) ,fTrkEovPBef(0) ,fTrkEovPAft(0) ,fdEdxBef(0) ,fdEdxAft(0) ,fIncpT(0) ,fIncpTM20(0) ,fInvmassLS(0) ,fInvmassULS(0) ,fInvmassLSmc(0) ,fInvmassULSmc(0) ,fInvmassLSreco(0) ,fInvmassULSreco(0) ,fInvmassLSmc0(0) ,fInvmassLSmc1(0) ,fInvmassLSmc2(0) ,fInvmassLSmc3(0) ,fInvmassULSmc0(0) ,fInvmassULSmc1(0) ,fInvmassULSmc2(0) ,fInvmassULSmc3(0) ,fOpeningAngleLS(0) ,fOpeningAngleULS(0) ,fPhotoElecPt(0) ,fPhoElecPt(0) ,fPhoElecPtM20(0) ,fPhoElecPtM20Mass(0) ,fSameElecPt(0) ,fSameElecPtM20(0) ,fSameElecPtM20Mass(0) ,fSemiElecPtM20(0) ,fTrackPtBefTrkCuts(0) ,fTrackPtAftTrkCuts(0) ,fTPCnsigma(0) ,fCent(0) ,fEleInfo(0) ,fElenSigma(0) /* ,fClsEBftTrigCut(0) ,fClsEAftTrigCut(0) ,fClsEAftTrigCut1(0) ,fClsEAftTrigCut2(0) ,fClsEAftTrigCut3(0) ,fClsEAftTrigCut4(0) ,fClsETime(0) ,fClsETime1(0) ,fTrigTimes(0) ,fCellCheck(0) */ ,fInputHFEMC(0) ,fInputAlle(0) ,fIncpTMChfe(0) ,fIncpTMChfeAll(0) ,fIncpTMCM20hfe(0) ,fIncpTMCM20hfeAll(0) ,fIncpTMCM20hfeCheck(0) ,fInputHFEMC_weight(0) ,fIncpTMCM20hfeCheck_weight(0) ,fIncpTMCpho(0) ,fIncpTMCM20pho(0) ,fPhoElecPtMC(0) ,fPhoElecPtMCM20(0) ,fPhoElecPtMCM20Mass(0) ,fSameElecPtMC(0) ,fSameElecPtMCM20(0) ,fSameElecPtMCM20Mass(0) ,fIncpTMCM20pho_pi0e(0) ,fPhoElecPtMCM20_pi0e(0) ,fSameElecPtMCM20_pi0e(0) ,fIncpTMCM20pho_eta(0) ,fPhoElecPtMCM20_eta(0) ,fSameElecPtMCM20_eta(0) ,fIncpTMCpho_pi0e_TPC(0) ,fPhoElecPtMC_pi0e_TPC(0) ,fSameElecPtMC_pi0e_TPC(0) ,fIncpTMCpho_eta_TPC(0) ,fPhoElecPtMC_eta_TPC(0) ,fSameElecPtMC_eta_TPC(0) ,CheckNclust(0) ,CheckNits(0) ,CheckDCA(0) ,Hpi0pTcheck(0) ,HETApTcheck(0) ,HphopTcheck(0) ,HDpTcheck(0) ,HBpTcheck(0) ,fpTCheck(0) ,fMomDtoE(0) ,fLabelCheck(0) ,fgeoFake(0) ,fFakeTrk0(0) ,fFakeTrk1(0) ,ftimingEle(0) ,fIncMaxE(0) ,fIncReco(0) ,fPhoReco(0) ,fSamReco(0) ,fIncRecoMaxE(0) ,fPhoRecoMaxE(0) ,fSamRecoMaxE(0) ,fPhoVertexReco_TPC(0) ,fPhoVertexReco_TPC_Invmass(0) ,fPhoVertexReco_EMCal(0) ,fPhoVertexReco_Invmass(0) ,fPhoVertexReco_step0(0) ,fPhoVertexReco_step1(0) ,fMatchV0_0(0) ,fMatchV0_1(0) ,fMatchMC_0(0) ,fMatchMC_1(0) ,fpair(0) ,fFakeRejection0(0) ,fFakeRejection1(0) ,fFakeRejection2(0) ,EopFake(0) ,EopTrue(0) ,MatchFake(0) ,MatchTrue(0) //,fnSigEtaCorr(NULL) { //Named constructor fPID = new AliHFEpid("hfePid"); fTrackCuts = new AliESDtrackCuts(); for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0; // Define input and output slots here // Input slot #0 works with a TChain DefineInput(0, TChain::Class()); // Output slot #0 id reserved by the base class for AOD // Output slot #1 writes into a TH1 container // DefineOutput(1, TH1I::Class()); DefineOutput(1, TList::Class()); // DefineOutput(3, TTree::Class()); } //________________________________________________________________________ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal") ,fESD(0) ,fMC(0) ,stack(0) ,fGeom(0) ,fOutputList(0) ,fqahist(0) ,fTrackCuts(0) ,fCuts(0) ,fIdentifiedAsOutInz(kFALSE) ,fPassTheEventCut(kFALSE) ,fRejectKinkMother(kFALSE) ,fmcData(kFALSE) ,fVz(0.0) ,fCFM(0) ,fPID(0) ,fPIDqa(0) ,fOpeningAngleCut(0.1) ,fMimpTassCut(0.5) ,fMimNsigassCut(-3) ,fInvmassCut(0) // no mass ,fSetMassConstraint(kTRUE) ,fSetMassWidthCut(kFALSE) ,fSetMassNonlinear(kFALSE) ,fSetKFpart(kTRUE) ,fNoEvents(0) ,fEMCAccE(0) ,hEMCAccE(0) ,fTrkpt(0) ,fTrkEovPBef(0) ,fTrkEovPAft(0) ,fdEdxBef(0) ,fdEdxAft(0) ,fIncpT(0) ,fIncpTM20(0) ,fInvmassLS(0) ,fInvmassULS(0) ,fInvmassLSmc(0) ,fInvmassULSmc(0) ,fInvmassLSreco(0) ,fInvmassULSreco(0) ,fInvmassLSmc0(0) ,fInvmassLSmc1(0) ,fInvmassLSmc2(0) ,fInvmassLSmc3(0) ,fInvmassULSmc0(0) ,fInvmassULSmc1(0) ,fInvmassULSmc2(0) ,fInvmassULSmc3(0) ,fOpeningAngleLS(0) ,fOpeningAngleULS(0) ,fPhotoElecPt(0) ,fPhoElecPt(0) ,fPhoElecPtM20(0) ,fPhoElecPtM20Mass(0) ,fSameElecPt(0) ,fSameElecPtM20(0) ,fSameElecPtM20Mass(0) ,fSemiElecPtM20(0) ,fTrackPtBefTrkCuts(0) ,fTrackPtAftTrkCuts(0) ,fTPCnsigma(0) ,fCent(0) ,fEleInfo(0) ,fElenSigma(0) /* ,fClsEBftTrigCut(0) ,fClsEAftTrigCut(0) ,fClsEAftTrigCut1(0) ,fClsEAftTrigCut2(0) ,fClsEAftTrigCut3(0) ,fClsEAftTrigCut4(0) ,fClsETime(0) ,fClsETime1(0) ,fTrigTimes(0) ,fCellCheck(0) */ ,fInputHFEMC(0) ,fInputAlle(0) ,fIncpTMChfe(0) ,fIncpTMChfeAll(0) ,fIncpTMCM20hfe(0) ,fIncpTMCM20hfeAll(0) ,fIncpTMCM20hfeCheck(0) ,fInputHFEMC_weight(0) ,fIncpTMCM20hfeCheck_weight(0) ,fIncpTMCpho(0) ,fIncpTMCM20pho(0) ,fPhoElecPtMC(0) ,fPhoElecPtMCM20(0) ,fPhoElecPtMCM20Mass(0) ,fSameElecPtMC(0) ,fSameElecPtMCM20(0) ,fSameElecPtMCM20Mass(0) ,fIncpTMCM20pho_pi0e(0) ,fPhoElecPtMCM20_pi0e(0) ,fSameElecPtMCM20_pi0e(0) ,fIncpTMCM20pho_eta(0) ,fPhoElecPtMCM20_eta(0) ,fSameElecPtMCM20_eta(0) ,fIncpTMCpho_pi0e_TPC(0) ,fPhoElecPtMC_pi0e_TPC(0) ,fSameElecPtMC_pi0e_TPC(0) ,fIncpTMCpho_eta_TPC(0) ,fPhoElecPtMC_eta_TPC(0) ,fSameElecPtMC_eta_TPC(0) ,CheckNclust(0) ,CheckNits(0) ,CheckDCA(0) ,Hpi0pTcheck(0) ,HETApTcheck(0) ,HphopTcheck(0) ,HDpTcheck(0) ,HBpTcheck(0) ,fpTCheck(0) ,fMomDtoE(0) ,fLabelCheck(0) ,fgeoFake(0) ,fFakeTrk0(0) ,fFakeTrk1(0) ,ftimingEle(0) ,fIncMaxE(0) ,fIncReco(0) ,fPhoReco(0) ,fSamReco(0) ,fIncRecoMaxE(0) ,fPhoRecoMaxE(0) ,fSamRecoMaxE(0) ,fPhoVertexReco_TPC(0) ,fPhoVertexReco_TPC_Invmass(0) ,fPhoVertexReco_EMCal(0) ,fPhoVertexReco_Invmass(0) ,fPhoVertexReco_step0(0) ,fPhoVertexReco_step1(0) ,fMatchV0_0(0) ,fMatchV0_1(0) ,fMatchMC_0(0) ,fMatchMC_1(0) ,fpair(0) ,fFakeRejection0(0) ,fFakeRejection1(0) ,fFakeRejection2(0) ,EopFake(0) ,EopTrue(0) ,MatchFake(0) ,MatchTrue(0) //,fnSigEtaCorr(NULL) { //Default constructor fPID = new AliHFEpid("hfePid"); fTrackCuts = new AliESDtrackCuts(); for(int i=0; i<7; i++)fnSigEtaCorr[i] = 0; // Constructor // Define input and output slots here // Input slot #0 works with a TChain DefineInput(0, TChain::Class()); // Output slot #0 id reserved by the base class for AOD // Output slot #1 writes into a TH1 container // DefineOutput(1, TH1I::Class()); DefineOutput(1, TList::Class()); //DefineOutput(3, TTree::Class()); } //_________________________________________ AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal() { //Destructor delete fOutputList; delete fGeom; delete fPID; delete fCuts; delete fCFM; delete fPIDqa; delete fTrackCuts; } //_________________________________________ void AliAnalysisTaskHFECal::UserExec(Option_t*) { //Main loop //Called for each event // create pointer to event fESD = dynamic_cast(InputEvent()); if (!fESD) { printf("ERROR: fESD not available\n"); return; } if(!fCuts){ AliError("HFE cuts not available"); return; } if(!fPID->IsInitialized()){ // Initialize PID with the given run number AliWarning("PID not initialised, get from Run no"); fPID->InitializePID(fESD->GetRunNumber()); } if(fmcData)fMC = MCEvent(); //AliStack* stack = NULL; if(fmcData && fMC)stack = fMC->Stack(); Float_t cent = -1.; AliCentrality *centrality = fESD->GetCentrality(); cent = centrality->GetCentralityPercentile("V0M"); //---- fill MC track info if(fmcData && fMC) { Int_t nParticles = stack->GetNtrack(); for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) { TParticle* particle = stack->Particle(iParticle); int fPDG = particle->GetPdgCode(); double mcZvertex = fMC->GetPrimaryVertex()->GetZ(); double pTMC = particle->Pt(); double proR = particle->R(); double etaMC = particle->Eta(); if(fabs(etaMC)>0.6)continue; Bool_t mcInDtoE= kFALSE; Bool_t mcInBtoE= kFALSE; Bool_t MChijing = fMC->IsFromBGEvent(iParticle); //if(!MChijing)printf("not MC hijing"); int iHijing = 1; if(!MChijing)iHijing = 0; double mcphoinfo[5]; mcphoinfo[0] = cent; mcphoinfo[1] = pTMC; mcphoinfo[2] = iHijing; if(fPDG==111)Hpi0pTcheck->Fill(mcphoinfo); if(fPDG==221)HETApTcheck->Fill(mcphoinfo); if(fabs(fPDG)==411 || fabs(fPDG)==413 || fabs(fPDG)==421 || fabs(fPDG)==423 || fabs(fPDG)==431)HDpTcheck->Fill(pTMC,iHijing); if(fabs(fPDG)==511 || fabs(fPDG)==513 || fabs(fPDG)==521 || fabs(fPDG)==523 || fabs(fPDG)==531)HBpTcheck->Fill(pTMC,iHijing); if(particle->GetFirstMother()>-1 && fPDG==22) { int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode(); if(parentPID==111 || parentPID==221)HphopTcheck->Fill(pTMC,iHijing); // pi0->g & eta->g } if(particle->GetFirstMother()>-1 && fabs(fPDG)==11) { int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode(); double pTMCparent = stack->Particle(particle->GetFirstMother())->Pt(); if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(fPDG)==11)mcInDtoE = kTRUE; if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(fPDG)==11)mcInBtoE = kTRUE; if((mcInBtoE || mcInDtoE) && fabs(mcZvertex)<10.0) { fInputHFEMC->Fill(cent,pTMC); double mcinfo[5]; mcinfo[0] = cent; mcinfo[1] = pTMC; mcinfo[2] = 0.0; mcinfo[3] = iHijing; mcinfo[4] = pTMCparent; fInputHFEMC_weight->Fill(mcinfo); } } if(proR<7.0 && fabs(fPDG)==11)fInputAlle->Fill(cent,pTMC); } } fNoEvents->Fill(0); Int_t fNOtrks = fESD->GetNumberOfTracks(); const AliESDVertex *pVtx = fESD->GetPrimaryVertex(); Double_t pVtxZ = -999; pVtxZ = pVtx->GetZ(); if(TMath::Abs(pVtxZ)>10) return; fNoEvents->Fill(1); if(fNOtrks<1) return; fNoEvents->Fill(2); AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse(); if(!pidResponse){ AliDebug(1, "Using default PID Response"); pidResponse = AliHFEtools::GetDefaultPID(kFALSE, fInputEvent->IsA() == AliAODEvent::Class()); } fPID->SetPIDResponse(pidResponse); fCFM->SetRecEventInfo(fESD); fCent->Fill(cent); //if(cent>90.) return; // Calorimeter info. FindTriggerClusters(); // make EMCAL array double maxE = 0.0; for(Int_t iCluster=0; iClusterGetNumberOfCaloClusters(); iCluster++) { AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster); if(clust && clust->IsEMCAL()) { double clustE = clust->E(); float emcx[3]; // cluster pos clust->GetPosition(emcx); TVector3 clustpos(emcx[0],emcx[1],emcx[2]); double emcphi = clustpos.Phi(); double emceta = clustpos.Eta(); double calInfo[5]; calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2(); //fEMCAccE->Fill(calInfo); hEMCAccE->Fill(cent,clustE); if(clustE>maxE)maxE = clustE; } } // Track loop for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { AliESDtrack* track = fESD->GetTrack(iTracks); if (!track) { printf("ERROR: Could not receive track %d\n", iTracks); continue; } //--- Get MC informtion int parentlabel = 99999; int parentPID = 99999; int grand_parentlabel = 99999; int grand_parentPID = 99999; Bool_t mcPho = kFALSE; Bool_t mcDtoE= kFALSE; Bool_t mcBtoE= kFALSE; Bool_t mcOrgPi0 = kFALSE; Bool_t mcOrgEta = kFALSE; double mcele = -1.; double mcpT = 0.0; double mcMompT = 0.0; //double mcGrandMompT = 0.0; double mcWeight = -10.0; double conv_proR = -1.0; int iHijing = 1; int mcLabel = -1; if(fmcData && fMC && stack) { Int_t label = TMath::Abs(track->GetLabel()); //mcLabel = track->GetLabel(); mcLabel = abs(track->GetLabel()); // check for conv. issue if(mcLabel>-1) { Bool_t MChijing = fMC->IsFromBGEvent(label); if(!MChijing)iHijing = 0; TParticle* particle = stack->Particle(label); int mcpid = particle->GetPdgCode(); mcpT = particle->Pt(); conv_proR = particle->R(); //printf("MCpid = %d",mcpid); if(particle->GetFirstMother()>-1) { //int parentlabel = particle->GetFirstMother(); //int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode(); //mcMompT = stack->Particle(particle->GetFirstMother())->Pt(); FindMother(particle, parentlabel, parentPID); mcMompT = stack->Particle(parentlabel)->Pt(); if((parentPID==22 || parentPID==111 || parentPID==221)&& fabs(mcpid)==11)mcPho = kTRUE; if((fabs(parentPID)==411 || fabs(parentPID)==413 || fabs(parentPID)==421 || fabs(parentPID)==423 || fabs(parentPID)==431)&& fabs(mcpid)==11)mcDtoE = kTRUE; if((fabs(parentPID)==511 || fabs(parentPID)==513 || fabs(parentPID)==521 || fabs(parentPID)==523 || fabs(parentPID)==531)&& fabs(mcpid)==11)mcBtoE = kTRUE; // make D->e pT correlation if(mcDtoE)fMomDtoE->Fill(mcpT,mcMompT); //cout << "check PID = " << parentPID << endl; //cout << "check pho = " << mcPho << endl; //cout << "check D or B = " << mcDtoE << endl; // pi->e (Dalitz) if(parentPID==111 && fabs(mcpid)==11 && mcMompT>0.0) { //cout << "find pi0->e " << endl; mcOrgPi0 = kTRUE; mcWeight = GetMCweight(mcMompT); } // eta->e (Dalitz) if(parentPID==221 && fabs(mcpid)==11 && mcMompT>0.0) { //cout << "find Eta->e " << endl; mcOrgEta = kTRUE; mcWeight = GetMCweightEta(mcMompT); } // access grand parent TParticle* particle_parent = stack->Particle(parentlabel); // get parent pointer //if(particle_parent->GetFirstMother()>-1 && parentPID==22 && fabs(mcpid)==11) // get grand parent g->e if(particle_parent->GetFirstMother()>-1 && (parentPID==22 || parentPID==111) && fabs(mcpid)==11) // get grand parent g->e { //int grand_parentPID = stack->Particle(particle_parent->GetFirstMother())->GetPdgCode(); //double pTtmp = stack->Particle(particle_parent->GetFirstMother())->Pt(); FindMother(particle_parent, grand_parentlabel, grand_parentPID); double mcGrandpT = stack->Particle(grand_parentlabel)->Pt(); if(grand_parentPID==111 && mcGrandpT>0.0) { // check eta->pi0 decay ! int grand2_parentlabel = 99999; int grand2_parentPID = 99999; TParticle* particle_grand = stack->Particle(grand_parentlabel); // get parent pointer FindMother(particle_grand, grand2_parentlabel, grand2_parentPID); if(grand2_parentPID==221) { //cout << "find Eta->e " << endl; double mcGrandpT2 = stack->Particle(grand2_parentlabel)->Pt(); mcOrgEta = kTRUE; mcWeight = GetMCweight(mcGrandpT2); mcMompT = mcGrandpT2; } else { //cout << "find pi0->e " << endl; mcOrgPi0 = kTRUE; mcWeight = GetMCweight(mcGrandpT); mcMompT = mcGrandpT; } } if(grand_parentPID==221 && mcGrandpT>0.0) { //cout << "find Eta->e " << endl; mcOrgEta = kTRUE; mcOrgPi0 = kFALSE; mcWeight = GetMCweightEta(mcGrandpT); mcMompT = mcGrandpT; } } } //cout << "===================="<e: " << mcele << endl; if(fabs(mcpid)==11 && mcBtoE)mcele= 2.; //cout << "check B->e: " << mcele << endl; if(fabs(mcpid)==11 && mcPho)mcele= 3.; //cout << "check Pho->e: " << mcele << endl; //cout << "check PID " << endl; if(fabs(mcpid)!=11) { //cout << "!= 11" << endl; //cout << mcpid << endl; } if(mcele==-1) { //cout << "mcele==-1" << endl; //cout << mcele << endl; //cout << mcpid << endl; } } // end of mcLabel>-1 } // <------ end of MC info. //cout << "Pi0 = " << mcOrgPi0 << " ; Eta = " << mcOrgEta << endl; //printf("weight = %f\n",mcWeight); if(TMath::Abs(track->Eta())>0.6) continue; //if(TMath::Abs(track->Pt()<2.5)) continue; if(TMath::Abs(track->Pt()<0.1)) continue; int nITS = track->GetNcls(0); fTrackPtBefTrkCuts->Fill(track->Pt()); /* UChar_t itsPixel = track->GetITSClusterMap(); cout << "nITS = " << nITS << endl; if(itsPixel & BIT(0))cout << "1st layer hit" << endl; if(itsPixel & BIT(1))cout << "2nd layer hit" << endl; if(itsPixel & BIT(2))cout << "3rd layer hit" << endl; if(itsPixel & BIT(3))cout << "4th layer hit" << endl; if(itsPixel & BIT(4))cout << "5th layer hit" << endl; */ // RecKine: ITSTPC cuts if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue; //RecKink if(fRejectKinkMother) { // Quick and dirty fix to reject both kink mothers and daughters if(track->GetKinkIndex(0) != 0) continue; } // RecPrim if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) continue; // HFEcuts: ITS layers cuts if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) continue; // HFE cuts: TPC PID cleanup if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue; if(mcPho && iHijing==0)fPhoVertexReco_step0->Fill(track->Pt(),conv_proR); // check MC vertex if(mcPho && iHijing==1)fPhoVertexReco_step1->Fill(track->Pt(),conv_proR); // check MC vertex int nTPCcl = track->GetTPCNcls(); //int nTPCclF = track->GetTPCNclsF(); // warnings //int nITS = track->GetNcls(0); fTrackPtAftTrkCuts->Fill(track->Pt()); Double_t mom = -999., eop=-999., pt = -999., dEdx=-999., fTPCnSigma=-10, phi=-999., eta=-999.; pt = track->Pt(); if(pt<2.5)continue; if(pt<0.1)continue; //Int_t charge = track->Charge(); fTrkpt->Fill(pt); mom = track->P(); phi = track->Phi(); eta = track->Eta(); float dca_xy; float dca_z; track->GetImpactParameters(dca_xy,dca_z); dEdx = track->GetTPCsignal(); fTPCnSigma = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(track, AliPID::kElectron) : 1000; //cout << "nSigma correctoon-----" << endl; //cout << "org = " << fTPCnSigma << endl; /* if(!fmcData) // nsigma eta correction { double nSigexpCorr = NsigmaCorrection(eta,cent); fTPCnSigma -= nSigexpCorr; } */ //cout << "correction = " << fTPCnSigma << endl; //--- track cluster match double ncells = -1.0; double m20 = -1.0; double m02 = -1.0; double disp = -1.0; double rmatch = -1.0; double nmatch = -1.0; //double oppstatus = 0.0; double emctof = 0.0; Bool_t MaxEmatch = kFALSE; Int_t clsId = track->GetEMCALcluster(); if (clsId>0){ AliESDCaloCluster *clust = fESD->GetCaloCluster(clsId); if(clust && clust->IsEMCAL()){ double clustE = clust->E(); if(clustE==maxE)MaxEmatch = kTRUE; eop = clustE/fabs(mom); //double clustT = clust->GetTOF(); ncells = clust->GetNCells(); m02 = clust->GetM02(); m20 = clust->GetM20(); disp = clust->GetDispersion(); double delphi = clust->GetTrackDx(); double deleta = clust->GetTrackDz(); rmatch = sqrt(pow(delphi,2)+pow(deleta,2)); nmatch = clust->GetNTracksMatched(); emctof = clust->GetTOF(); //cout << "emctof = " << emctof << endl; cout << "eop org = "<< eop << endl; double eoporg = eop; if(fmcData) { double mceopcorr = MCEopMeanCorrection(pt,cent); eop += mceopcorr; } cout << "eop corr = " << eop << endl; double valdedx[16]; valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma; valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = dca_xy, valdedx[8] = dca_z; valdedx[9] = m20; valdedx[10] = mcpT; valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = eoporg; valdedx[14] = nTPCcl; valdedx[15] = mcele; fEleInfo->Fill(valdedx); } } //Get Cal info PID response /* double eop2; double ss[4]; Double_t nSigmaEop = fPID->GetPIDResponse()->NumberOfSigmasEMCAL(track,AliPID::kElectron,eop2,ss); if(fTPCnSigma>-1.5 && fTPCnSigma<3.0 && nITS>2.5 && nTPCcl>100) { double valEop[3]; valEop[0] = cent; valEop[1] = pt; valEop[2] = nSigmaEop; fElenSigma->Fill(valEop); } */ // --- tarck cut & e ID if(nITS<2.5)continue; if(nTPCcl<100)continue; CheckNclust->Fill(nTPCcl); CheckNits->Fill(nITS); CheckDCA->Fill(dca_xy,dca_z); // check production vertex of photons fdEdxBef->Fill(mom,fTPCnSigma); fTPCnsigma->Fill(mom,fTPCnSigma); if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop); //--- track accepted by HFE Int_t pidpassed = 1; AliHFEpidObject hfetrack; hfetrack.SetAnalysisType(AliHFEpidObject::kESDanalysis); hfetrack.SetRecTrack(track); double binct = 10.5; if((0.0< cent) && (cent<5.0)) binct = 0.5; if((5.0< cent) && (cent<10.0)) binct = 1.5; if((10.0< cent) && (cent<20.0)) binct = 2.5; if((20.0< cent) && (cent<30.0)) binct = 3.5; if((30.0< cent) && (cent<40.0)) binct = 4.5; if((40.0< cent) && (cent<50.0)) binct = 5.5; if((50.0< cent) && (cent<60.0)) binct = 6.5; if((60.0< cent) && (cent<70.0)) binct = 7.5; if((70.0< cent) && (cent<80.0)) binct = 8.5; if((80.0< cent) && (cent<90.0)) binct = 9.5; if((90.0< cent) && (cent<100.0)) binct = 10.5; hfetrack.SetCentrality((int)binct); //added hfetrack.SetPbPb(); if(!fPID->IsSelected(&hfetrack, NULL, "", fPIDqa)) pidpassed = 0; if(pidpassed==0) continue; // nSigma rejection //--- photonic ID Bool_t fFlagPhotonicTPC = kFALSE; Bool_t fFlagConvinatTPC = kFALSE; SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicTPC,fFlagConvinatTPC,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta,0); //--- check reco eff. with TPC double phoval[5]; phoval[0] = cent; phoval[1] = pt; phoval[2] = fTPCnSigma; phoval[3] = iHijing; phoval[4] = mcMompT; if((fTPCnSigma >= -5.0 && fTPCnSigma <= 5) && (mcOrgPi0 || mcOrgEta)) { if(iHijing==1)mcWeight = 1.0; if(mcOrgPi0) { fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight); if(fFlagPhotonicTPC) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight); if(fFlagConvinatTPC) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight); fPhoVertexReco_TPC->Fill(track->Pt(),conv_proR,mcWeight); // check MC vertex if(fFlagPhotonicTPC)fPhoVertexReco_TPC_Invmass->Fill(track->Pt(),conv_proR,mcWeight); // check MC vertex } if(mcOrgEta) { fIncpTMCpho_eta_TPC->Fill(phoval,mcWeight); if(fFlagPhotonicTPC) fPhoElecPtMC_eta_TPC->Fill(phoval,mcWeight); if(fFlagConvinatTPC) fSameElecPtMC_eta_TPC->Fill(phoval,mcWeight); } } //--- matching check double emcphimim = 1.396; double emcphimax = 3.14; if(phi>emcphimim && phiFill(pt); // data if(mcele>2.1)fMatchMC_0->Fill(pt,mcWeight); // MC if(eop>=0.0) // have a match { if(fFlagPhotonicTPC)fMatchV0_1->Fill(pt); if(mcele>2.1)fMatchMC_1->Fill(pt,mcWeight); } } // check fake rejection if(mcOrgPi0 || mcOrgEta) { double phiacc0 = 80.0/180.0*acos(-1); double phiacc1 = 180.0/180.0*acos(-1); int TrStat = 0; if(phi>phiacc0 && phiGetLabel()>0)TrStat = 1; fFakeRejection0->Fill(TrStat,pt,mcWeight); if(eop>-1.0)fFakeRejection1->Fill(TrStat,pt,mcWeight); // have match if(eop>0.9 && eop<1.3)fFakeRejection2->Fill(TrStat,pt,mcWeight); // have PID if(TrStat==0) { EopFake->Fill(pt,eop,mcWeight); MatchFake->Fill(pt,rmatch,mcWeight); } else { EopTrue->Fill(pt,eop,mcWeight); MatchTrue->Fill(pt,rmatch,mcWeight); } } } //+++++++ E/p cut ++++++++++++++++ if(eop<0.9 || eop>1.3)continue; Bool_t fFlagPhotonicElec = kFALSE; Bool_t fFlagConvinatElec = kFALSE; SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta,1); fTrkEovPAft->Fill(pt,eop); fdEdxAft->Fill(mom,fTPCnSigma); // Fill real data fIncpT->Fill(cent,pt); if(fFlagPhotonicElec) fPhoElecPt->Fill(cent,pt); if(fFlagConvinatElec) fSameElecPt->Fill(cent,pt); if(m20>0.0 && m20<0.3) { fIncpTM20->Fill(cent,pt); ftimingEle->Fill(pt,emctof); if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt); if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt); if(!fFlagPhotonicElec) fSemiElecPtM20->Fill(cent,pt); } // MC // check label for electron candidiates int idlabel = 1; if(mcLabel==0)idlabel = 0; fLabelCheck->Fill(pt,idlabel); if(mcLabel==0)fgeoFake->Fill(phi,eta); if(mcLabel<0 && m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3) { fFakeTrk0->Fill(cent,pt); } if(mcele>-1) // select MC electrons { fIncpTMChfeAll->Fill(cent,pt); if(m20>0.0 && m20<0.3)fIncpTMCM20hfeAll->Fill(cent,pt); if(m20>0.0 && m20<0.3 && fTPCnSigma>-1 && fTPCnSigma<3)fFakeTrk1->Fill(cent,pt); if(mcBtoE || mcDtoE) // select B->e & D->e { fIncpTMChfe->Fill(cent,pt); if(m20>0.0 && m20<0.3) { //cout << "MC label = " << mcLabel << endl; fIncpTMCM20hfe->Fill(cent,pt); fIncpTMCM20hfeCheck->Fill(cent,mcpT); fIncpTMCM20hfeCheck_weight->Fill(phoval); } } if(mcPho) // select photonic electrons { fIncpTMCpho->Fill(phoval); if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval); if(fFlagConvinatElec) fSameElecPtMC->Fill(phoval); if(m20>0.0 && m20<0.3) { fIncpTMCM20pho->Fill(phoval); if(fFlagPhotonicElec) fPhoElecPtMCM20->Fill(phoval); if(fFlagConvinatElec) fSameElecPtMCM20->Fill(phoval); // pi0->g->e if(mcWeight>-1) { if(iHijing==1)mcWeight = 1.0; if(mcOrgPi0) { fIncpTMCM20pho_pi0e->Fill(phoval,mcWeight); if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval,mcWeight); if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval,mcWeight); // check production vertex //fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR); //if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR); fPhoVertexReco_EMCal->Fill(track->Pt(),conv_proR,mcWeight); if(fFlagPhotonicElec)fPhoVertexReco_Invmass->Fill(track->Pt(),conv_proR,mcWeight); } // --- eta if(mcOrgEta) { fIncpTMCM20pho_eta->Fill(phoval,mcWeight); if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight); if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight); } // check production vertex } } } } } PostData(1, fOutputList); } //_________________________________________ void AliAnalysisTaskHFECal::UserCreateOutputObjects() { //--- Check MC //Bool_t mcData = kFALSE; if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) { fmcData = kTRUE; printf("+++++ MC Data available"); } if(fmcData) { printf("++++++++= MC analysis \n"); } else { printf("++++++++ real data analysis \n"); } printf("+++++++ QA hist %d",fqahist); //---- Geometry fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1"); //--------Initialize PID //fPID->SetHasMCData(kFALSE); fPID->SetHasMCData(fmcData); if(!fPID->GetNumberOfPIDdetectors()) { fPID->AddDetector("TPC", 0); //fPID->AddDetector("EMCAL", 1); <--- apply PID selection in task } Double_t params[4]; const char *cutmodel; cutmodel = "pol0"; params[0] = -1.0; //sigma min double maxnSig = 3.0; if(fmcData) { params[0] = -5.0; //sigma min maxnSig = 5.0; } for(Int_t a=0;a<11;a++)fPID->ConfigureTPCcentralityCut(a,cutmodel,params,maxnSig); fPID->SortDetectors(); fPIDqa = new AliHFEpidQAmanager(); fPIDqa->Initialize(fPID); //------- fcut -------------- fCuts = new AliHFEcuts(); fCuts->CreateStandardCuts(); //fCuts->SetMinNClustersTPC(100); fCuts->SetMinNClustersTPC(90); fCuts->SetMinRatioTPCclusters(0.6); fCuts->SetTPCmodes(AliHFEextraCuts::kFound, AliHFEextraCuts::kFoundOverFindable); //fCuts->SetMinNClustersITS(3); fCuts->SetMinNClustersITS(2); fCuts->SetProductionVertex(0,50,0,50); fCuts->SetCutITSpixel(AliHFEextraCuts::kAny); fCuts->SetCheckITSLayerStatus(kFALSE); fCuts->SetVertexRange(10.); fCuts->SetTOFPIDStep(kFALSE); //fCuts->SetPtRange(2, 50); fCuts->SetPtRange(0.1, 50); //fCuts->SetMaxImpactParam(3.,3.); fCuts->SetMaxImpactParam(2.4,3.2); // standard in 2011 //--------Initialize correction Framework and Cuts fCFM = new AliCFManager; const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack; fCFM->SetNStepParticle(kNcutSteps); for(Int_t istep = 0; istep < kNcutSteps; istep++) fCFM->SetParticleCutsList(istep, NULL); if(!fCuts){ AliWarning("Cuts not available. Default cuts will be used"); fCuts = new AliHFEcuts; fCuts->CreateStandardCuts(); } fCuts->Initialize(fCFM); //---------Output Tlist fOutputList = new TList(); fOutputList->SetOwner(); fOutputList->Add(fPIDqa->MakeList("PIDQA")); fNoEvents = new TH1F("fNoEvents","",4,-0.5,3.5) ; fOutputList->Add(fNoEvents); Int_t binsE[5] = {250, 100, 1000, 200, 10}; Double_t xminE[5] = {1.0, -1, 0.0, 0, -0.5}; Double_t xmaxE[5] = {3.5, 1, 100.0, 100, 9.5}; fEMCAccE = new THnSparseD("fEMCAccE","EMC acceptance & E;#eta;#phi;Energy;Centrality;trugCondition;",5,binsE,xminE,xmaxE); //if(fqahist==1)fOutputList->Add(fEMCAccE); hEMCAccE = new TH2F("hEMCAccE","Cluster Energy",200,0,100,100,0,20); fOutputList->Add(hEMCAccE); fTrkpt = new TH1F("fTrkpt","track pt",100,0,50); fOutputList->Add(fTrkpt); fTrackPtBefTrkCuts = new TH1F("fTrackPtBefTrkCuts","track pt before track cuts",100,0,50); fOutputList->Add(fTrackPtBefTrkCuts); fTrackPtAftTrkCuts = new TH1F("fTrackPtAftTrkCuts","track pt after track cuts",100,0,50); fOutputList->Add(fTrackPtAftTrkCuts); fTPCnsigma = new TH2F("fTPCnsigma", "TPC - n sigma",100,0,50,200,-10,10); fOutputList->Add(fTPCnsigma); fTrkEovPBef = new TH2F("fTrkEovPBef","track E/p before HFE pid",100,0,50,100,0,2); fOutputList->Add(fTrkEovPBef); fTrkEovPAft = new TH2F("fTrkEovPAft","track E/p after HFE pid",100,0,50,100,0,2); fOutputList->Add(fTrkEovPAft); fdEdxBef = new TH2F("fdEdxBef","track dEdx vs p before HFE pid",100,0,50,200,-10,10); fOutputList->Add(fdEdxBef); fdEdxAft = new TH2F("fdEdxAft","track dEdx vs p after HFE pid",100,0,50,200,-10,10); fOutputList->Add(fdEdxAft); fIncpT = new TH2F("fIncpT","HFE pid electro vs. centrality",200,0,100,100,0,50); fOutputList->Add(fIncpT); fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50); fOutputList->Add(fIncpTM20); Int_t nBinspho[9] = { 10, 30, 600, 60, 50, 4, 40, 8, 30}; Double_t minpho[9] = { 0., 0., -0.1, 40, 0, -0.5, 0,-1.5, 0}; Double_t maxpho[9] = {100., 30., 0.5, 100, 1, 3.5, 2, 6.5, 30}; fInvmassLS = new THnSparseD("fInvmassLS", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho); if(fqahist==1)fOutputList->Add(fInvmassLS); fInvmassULS = new THnSparseD("fInvmassULS", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho); if(fqahist==1)fOutputList->Add(fInvmassULS); fInvmassLSmc = new THnSparseD("fInvmassLSmc", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; Mcele;", 9, nBinspho,minpho, maxpho); if(fqahist==1)fOutputList->Add(fInvmassLSmc); fInvmassULSmc = new THnSparseD("fInvmassULSmc", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2); nSigma; angle; m20cut; eop; MCele", 9, nBinspho,minpho, maxpho); if(fqahist==1)fOutputList->Add(fInvmassULSmc); fInvmassLSreco = new TH2D("fInvmassLSreco", "Inv mass of LS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassLSreco->Sumw2(); fOutputList->Add(fInvmassLSreco); fInvmassULSreco = new TH2D("fInvmassULSreco", "Inv mass of ULS (e,e) reco; cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassULSreco->Sumw2(); fOutputList->Add(fInvmassULSreco); fInvmassLSmc0 = new TH2D("fInvmassLSmc0", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassLSmc0->Sumw2(); fOutputList->Add(fInvmassLSmc0); fInvmassLSmc1 = new TH2D("fInvmassLSmc1", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassLSmc1->Sumw2(); fOutputList->Add(fInvmassLSmc1); fInvmassLSmc2 = new TH2D("fInvmassLSmc2", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassLSmc2->Sumw2(); fOutputList->Add(fInvmassLSmc2); fInvmassLSmc3 = new TH2D("fInvmassLSmc3", "Inv mass of LS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassLSmc3->Sumw2(); fOutputList->Add(fInvmassLSmc3); fInvmassULSmc0 = new TH2D("fInvmassULSmc0", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassULSmc0->Sumw2(); fOutputList->Add(fInvmassULSmc0); fInvmassULSmc1 = new TH2D("fInvmassULSmc1", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassULSmc1->Sumw2(); fOutputList->Add(fInvmassULSmc1); fInvmassULSmc2 = new TH2D("fInvmassULSmc2", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassULSmc2->Sumw2(); fOutputList->Add(fInvmassULSmc2); fInvmassULSmc3 = new TH2D("fInvmassULSmc3", "Inv mass of ULS (e,e); cent; p_{T} (GeV/c); mass(GeV/c^2)",20,0,20,600,-0.1,0.5 ); fInvmassULSmc3->Sumw2(); fOutputList->Add(fInvmassULSmc3); fOpeningAngleLS = new TH1F("fOpeningAngleLS","Opening angle for LS pairs",100,0,1); fOutputList->Add(fOpeningAngleLS); fOpeningAngleULS = new TH1F("fOpeningAngleULS","Opening angle for ULS pairs",100,0,1); fOutputList->Add(fOpeningAngleULS); fPhotoElecPt = new TH1F("fPhotoElecPt", "photonic electron pt",100,0,50); fOutputList->Add(fPhotoElecPt); fPhoElecPt = new TH2F("fPhoElecPt", "Pho-inclusive electron pt",200,0,100,100,0,50); fOutputList->Add(fPhoElecPt); fPhoElecPtM20 = new TH2F("fPhoElecPtM20", "Pho-inclusive electron pt with M20",200,0,100,100,0,50); fOutputList->Add(fPhoElecPtM20); fPhoElecPtM20Mass = new TH2F("fPhoElecPtM20Mass", "Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50); fPhoElecPtM20Mass->Sumw2(); fOutputList->Add(fPhoElecPtM20Mass); fSameElecPt = new TH2F("fSameElecPt", "Same-inclusive electron pt",200,0,100,100,0,50); fOutputList->Add(fSameElecPt); fSameElecPtM20 = new TH2F("fSameElecPtM20", "Same-inclusive electron pt with M20",200,0,100,100,0,50); fOutputList->Add(fSameElecPtM20); fSameElecPtM20Mass = new TH2F("fSameElecPtM20Mass", "Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50); fSameElecPtM20Mass->Sumw2(); fOutputList->Add(fSameElecPtM20Mass); fSemiElecPtM20 = new TH2F("fSemiElecPtM20", "Semi-inclusive electron pt with M20",200,0,100,100,0,50); fOutputList->Add(fSemiElecPtM20); fCent = new TH1F("fCent","Centrality",200,0,100) ; fOutputList->Add(fCent); // Make common binning const Double_t kMinP = 0.; const Double_t kMaxP = 20.; //+++ 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, cell, M02, M20, Disp, Centrality, select) // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, phi, eta, Sig, e/p, ,match, dca_xy, dca_z, M20, Disp, Centrality, select) Int_t nBins[16] = { 100, 7, 60, 20, 90, 100, 25, 60, 60, 100, 40, 10, 250, 100, 100, 8}; Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -5.0, 0, 0, -3.0, -3.0, 0.0, 0.0, 0, 0, 0, 80, -1.5}; Double_t max[16] = {kMaxP, 6.5, 4.0, 1.0, 4.0, 2.0, 0.05, 3.0, 3.0, 1.0, 20.0, 100, 100, 2.0, 180, 6.5}; //fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;clsF;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max); fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;DCA_xy;DCA_z;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max); //if(fqahist==1)fOutputList->Add(fEleInfo); // Make common binning Int_t nBinsEop[3] = { 10, 50, 100}; Double_t minEop[3] = { 0, 0, -5}; Double_t maxEop[3] = {100, 50, 5}; fElenSigma= new THnSparseD("fElenSigma", "Electron nSigma; cent; pT [GeV/c]; nSigma", 3, nBinsEop, minEop, maxEop); fOutputList->Add(fElenSigma); //<--- Trigger info /* fClsEBftTrigCut = new TH1F("fClsEBftTrigCut","cluster E before trigger selection",1000,0,100); fOutputList->Add(fClsEBftTrigCut); fClsEAftTrigCut = new TH1F("fClsEAftTrigCut","cluster E if cls has 0 trigcut channel",1000,0,100); fOutputList->Add(fClsEAftTrigCut); fClsEAftTrigCut1 = new TH1F("fClsEAftTrigCut1","cluster E if cls with trig channel",1000,0,100); fOutputList->Add(fClsEAftTrigCut1); fClsEAftTrigCut2 = new TH1F("fClsEAftTrigCut2","cluster E if cls with trigcut channel",1000,0,100); fOutputList->Add(fClsEAftTrigCut2); fClsEAftTrigCut3 = new TH1F("fClsEAftTrigCut3","cluster E if cls with trigcut channel + nCell>Ecorrect",1000,0,100); fOutputList->Add(fClsEAftTrigCut3); fClsEAftTrigCut4 = new TH1F("fClsEAftTrigCut4","cluster E if cls with trigcut channel + nCell>Ecorrect + cls time cut",1000,0,100); fOutputList->Add(fClsEAftTrigCut4); fClsETime = new TH2F("fClsETime", "Cls time vs E; E; time;",1000,0,100,1000,-0.0000002,0.0000009); fOutputList->Add(fClsETime); fClsETime1 = new TH2F("fClsETime1", "Cls time vs E if cls contains trigger channel; E; time;",1000,0,100,1000,-0.0000002,0.0000009); fOutputList->Add(fClsETime1); fTrigTimes = new TH1F("fTrigTimes", "Trigger time; time; N;",25,0,25); fOutputList->Add(fTrigTimes); fCellCheck = new TH2F("fCellCheck", "Cell vs E; E GeV; Cell ID",10,6,26,12000,0,12000); fOutputList->Add(fCellCheck); */ //<---------- MC fInputHFEMC = new TH2F("fInputHFEMC","Input MC HFE pid electro vs. centrality",200,0,100,100,0,50); fOutputList->Add(fInputHFEMC); fInputAlle = new TH2F("fInputAlle","Input MC electro vs. centrality",200,0,100,100,0,50); fOutputList->Add(fInputAlle); fIncpTMChfe = new TH2F("fIncpTMChfe","MC HFE pid electro vs. centrality",200,0,100,100,0,50); fOutputList->Add(fIncpTMChfe); fIncpTMChfeAll = new TH2F("fIncpTMChfeAll","MC Alle pid electro vs. centrality",200,0,100,100,0,50); fOutputList->Add(fIncpTMChfeAll); fIncpTMCM20hfe = new TH2F("fIncpTMCM20hfe","MC HFE pid electro vs. centrality with M20",200,0,100,100,0,50); fOutputList->Add(fIncpTMCM20hfe); fIncpTMCM20hfeAll = new TH2F("fIncpTMCM20hfeAll","MC Alle pid electro vs. centrality with M20",200,0,100,100,0,50); fOutputList->Add(fIncpTMCM20hfeAll); fIncpTMCM20hfeCheck = new TH2F("fIncpTMCM20hfeCheck","MC HFE pid electro vs. centrality with M20 Check",200,0,100,100,0,50); fOutputList->Add(fIncpTMCM20hfeCheck); Int_t nBinspho2[5] = { 200, 100, 7, 3, 700}; Double_t minpho2[5] = { 0., 0., -2.5, -0.5, 0.}; Double_t maxpho2[5] = {100., 50., 4.5, 2.5, 70.}; fInputHFEMC_weight = new THnSparseD("fInputHFEMC_weight", "MC HFE electron pt",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fInputHFEMC_weight); fIncpTMCM20hfeCheck_weight = new THnSparseD("fIncpTMCM20hfeCheck_weight", "HFE electron pt with M20",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fIncpTMCM20hfeCheck_weight); fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fIncpTMCpho); fIncpTMCM20pho = new THnSparseD("fIncpTMCM20pho","MC Pho pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fIncpTMCM20pho); fPhoElecPtMC = new THnSparseD("fPhoElecPtMC", "MC Pho-inclusive electron pt",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fPhoElecPtMC); fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fPhoElecPtMCM20); fPhoElecPtMCM20Mass = new TH2D("fPhoElecPtMCM20Mass", "MC Pho-inclusive electron pt with M20 Mass",200,0,100,100,0,50); fOutputList->Add(fPhoElecPtMCM20Mass); fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fSameElecPtMC); fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fSameElecPtMCM20); fSameElecPtMCM20Mass = new TH2D("fSameElecPtMCM20Mass", "MC Same-inclusive electron pt with M20 Mass",200,0,100,100,0,50); fOutputList->Add(fSameElecPtMCM20Mass); fIncpTMCM20pho_pi0e = new THnSparseD("fIncpTMCM20pho_pi0e","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2); fIncpTMCM20pho_pi0e->Sumw2(); fOutputList->Add(fIncpTMCM20pho_pi0e); fPhoElecPtMCM20_pi0e = new THnSparseD("fPhoElecPtMCM20_pi0e", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2); fPhoElecPtMCM20_pi0e->Sumw2(); fOutputList->Add(fPhoElecPtMCM20_pi0e); fSameElecPtMCM20_pi0e = new THnSparseD("fSameElecPtMCM20_pi0e", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2); fSameElecPtMCM20_pi0e->Sumw2(); fOutputList->Add(fSameElecPtMCM20_pi0e); // fIncpTMCM20pho_eta = new THnSparseD("fIncpTMCM20pho_eta","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2); fIncpTMCM20pho_eta->Sumw2(); fOutputList->Add(fIncpTMCM20pho_eta); fPhoElecPtMCM20_eta = new THnSparseD("fPhoElecPtMCM20_eta", "MC Pho-inclusive electron pt with M20 pi0->e",5,nBinspho2,minpho2,maxpho2); fPhoElecPtMCM20_eta->Sumw2(); fOutputList->Add(fPhoElecPtMCM20_eta); fSameElecPtMCM20_eta = new THnSparseD("fSameElecPtMCM20_eta", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2); fSameElecPtMCM20_eta->Sumw2(); fOutputList->Add(fSameElecPtMCM20_eta); // ------------ fIncpTMCpho_pi0e_TPC = new THnSparseD("fIncpTMCpho_pi0e_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2); fIncpTMCpho_pi0e_TPC->Sumw2(); fOutputList->Add(fIncpTMCpho_pi0e_TPC); fPhoElecPtMC_pi0e_TPC = new THnSparseD("fPhoElecPtMC_pi0e_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2); fPhoElecPtMC_pi0e_TPC->Sumw2(); fOutputList->Add(fPhoElecPtMC_pi0e_TPC); fSameElecPtMC_pi0e_TPC = new THnSparseD("fSameElecPtMC_pi0e_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2); fSameElecPtMC_pi0e_TPC->Sumw2(); fOutputList->Add(fSameElecPtMC_pi0e_TPC); fIncpTMCpho_eta_TPC = new THnSparseD("fIncpTMCpho_eta_TPC","MC Pho pi0->e pid electro vs. centrality with M20",5,nBinspho2,minpho2,maxpho2); fIncpTMCpho_eta_TPC->Sumw2(); fOutputList->Add(fIncpTMCpho_eta_TPC); fPhoElecPtMC_eta_TPC = new THnSparseD("fPhoElecPtMC_eta_TPC", "MC Pho-inclusive electron pt with pi0->e",5,nBinspho2,minpho2,maxpho2); fPhoElecPtMC_eta_TPC->Sumw2(); fOutputList->Add(fPhoElecPtMC_eta_TPC); fSameElecPtMC_eta_TPC = new THnSparseD("fSameElecPtMC_eta_TPC", "MC Same-inclusive electron pt pi0->e",5,nBinspho2,minpho2,maxpho2); fSameElecPtMC_eta_TPC->Sumw2(); fOutputList->Add(fSameElecPtMC_eta_TPC); //------------- CheckNclust = new TH1D("CheckNclust","cluster check",200,0,200); fOutputList->Add(CheckNclust); CheckNits = new TH1D("CheckNits","ITS cluster check",8,-0.5,7.5); fOutputList->Add(CheckNits); CheckDCA = new TH2D("CheckDCA","DCA check",200,-5,5,200,-5,5); fOutputList->Add(CheckDCA); /* Hpi0pTcheck = new TH2D("Hpi0pTcheck","Pi0 pT from Hijing",100,0,50,3,-0.5,2.5); fOutputList->Add(Hpi0pTcheck); HETApTcheck = new TH2D("HETApTcheck","Eta pT from Hijing",100,0,50,3,-0.5,2.5); fOutputList->Add(HETApTcheck); */ Int_t nBinspho3[3] = { 200, 100, 3}; Double_t minpho3[3] = { 0., 0., -0.5}; Double_t maxpho3[3] = {100., 50., 2.5}; Hpi0pTcheck = new THnSparseD("Hpi0pTcheck","Pi0 pT from Hijing",3,nBinspho3,minpho3,maxpho3); fOutputList->Add(Hpi0pTcheck); HETApTcheck = new THnSparseD("HETApTcheck","Eta pT from Hijing",3,nBinspho3,minpho3,maxpho3); fOutputList->Add(HETApTcheck); //-- HphopTcheck = new TH2D("HphopTcheck","Pho pT from Hijing",100,0,50,3,-0.5,2.5); fOutputList->Add(HphopTcheck); // HDpTcheck = new TH2D("HDpTcheck","D pT from Hijing",100,0,50,3,-0.5,2.5); fOutputList->Add(HDpTcheck); HBpTcheck = new TH2D("HBpTcheck","B pT from Hijing",100,0,50,3,-0.5,2.5); fOutputList->Add(HBpTcheck); // fpTCheck = new TH1D("fpTCheck","pT check",500,0,50); fOutputList->Add(fpTCheck); fMomDtoE = new TH2D("fMomDtoE","D->E pT correlations;e p_{T} GeV/c;D p_{T} GeV/c",400,0,40,400,0,40); fOutputList->Add(fMomDtoE); fLabelCheck = new TH2D("fLabelCheck","MC label",50,0,50,5,-1.5,3.5); fOutputList->Add(fLabelCheck); fgeoFake = new TH2D("fgeoFake","Label==0 eta and phi",628,0,6.28,200,-1,1); fOutputList->Add(fgeoFake); fFakeTrk0 = new TH2D("fFakeTrk0","fake trakcs",10,0,100,20,0,20); fOutputList->Add(fFakeTrk0); fFakeTrk1 = new TH2D("fFakeTrk1","true all e a.f. eID",10,0,100,20,0,20); fOutputList->Add(fFakeTrk1); ftimingEle = new TH2D("ftimingEle","electron TOF",100,0,20,100,1e-7,1e-6); fOutputList->Add(ftimingEle); // eta correction // note: parameters 01/31new.TPCnSigmaEtaDep // 70-90 delta_eta = 0.2 double etaval[12] = {-0.55,-0.45,-0.35,-0.25,-0.15,-0.05,0.05,0.15,0.25,0.35,0.45,0.55}; double corr0[12]= {-0.569177,-0.528844,-0.391979,-0.165494,0.0283495,0.156171,0.266353,0.13103,-0.0250842,-0.274089,-0.45481,-0.536291}; // 0-10 (done) double corr1[12]= {-0.404742,-0.278953,-0.218069,0.00139927,0.191412,0.354403,0.524594,0.341778,0.244199,-0.112146,-0.160692,-0.352832}; // 10-20 (done) double corr2[12] = {-0.306007,-0.16821,-0.0248635,0.202233,0.447051,0.497197,0.712251,0.433482,0.337907,0.168426,-0.0693229,-0.0728351}; // 20-30 (done) double corr3[12] = {-0.13884,-0.0503553,0.104403,0.389773,0.50697,0.539048,0.751642,0.655636,0.518563,0.308156,0.0361159,-0.0491439}; // 30-40 (done) double corr4[12] = {-0.0319431,0.0808711,0.208774,0.443217,0.557762,0.61453,0.889519,0.808282,0.620394,0.267092,0.15241,-0.0458664}; // 40-50 (done) double corr5[12] = {-0.130625,0.0189124,0.190344,0.467431,0.546353,0.672251,0.731541,0.802101,0.437108,0.294081,0.193682,0.159074}; // 50-70(done) double corr6[12] = {0.0600197,0.0600197,0.358366,0.358366,0.973734,0.973734,0.759812,0.759812,0.667861,0.667861,0.415635,0.415635}; // 70-90(done) fnSigEtaCorr[0] = new TGraphErrors(12,etaval,corr0); // 0-10 fnSigEtaCorr[1] = new TGraphErrors(12,etaval,corr1); // 10-20 fnSigEtaCorr[2] = new TGraphErrors(12,etaval,corr2); // 20-30 fnSigEtaCorr[3] = new TGraphErrors(12,etaval,corr3); // 30-40 fnSigEtaCorr[4] = new TGraphErrors(12,etaval,corr4); // 40-50 fnSigEtaCorr[5] = new TGraphErrors(12,etaval,corr5); // 50-70 fnSigEtaCorr[6] = new TGraphErrors(12,etaval,corr6); // 70-90 fIncMaxE = new TH2D("fIncMaxE","Inc",10,0,100,10,0,100); fOutputList->Add(fIncMaxE); fIncReco = new TH2D("fIncReco","Inc",10,0,100,100,0,500); fOutputList->Add(fIncReco); fPhoReco = new TH2D("fPhoReco","Pho",10,0,100,100,0,500); fOutputList->Add(fPhoReco); fSamReco = new TH2D("fSamReco","Same",10,0,100,100,0,500); fOutputList->Add(fSamReco); fIncRecoMaxE = new TH2D("fIncRecoMaxE","Inc",10,0,100,100,0,500); fOutputList->Add(fIncRecoMaxE); fPhoRecoMaxE = new TH2D("fPhoRecoMaxE","Pho",10,0,100,100,0,500); fOutputList->Add(fPhoRecoMaxE); fSamRecoMaxE = new TH2D("fSamRecoMaxE","Same",10,0,100,100,0,500); fOutputList->Add(fSamRecoMaxE); fPhoVertexReco_TPC = new TH2D("fPhoVertexReco_TPC","photon production Vertex mass selection TPC",40,0,20,250,0,50); fPhoVertexReco_TPC->Sumw2(); fOutputList->Add(fPhoVertexReco_TPC); fPhoVertexReco_TPC_Invmass = new TH2D("fPhoVertexReco_TPC_Invmass","photon production Vertex mass selection TPC Invmass",40,0,20,250,0,50); fPhoVertexReco_TPC_Invmass->Sumw2(); fOutputList->Add(fPhoVertexReco_TPC_Invmass); fPhoVertexReco_EMCal = new TH2D("fPhoVertexReco_EMCal","photon production Vertex mass selection",40,0,20,250,0,50); fPhoVertexReco_EMCal->Sumw2(); fOutputList->Add(fPhoVertexReco_EMCal); fPhoVertexReco_Invmass = new TH2D("fPhoVertexReco_Invmass","photon production Vertex mass selection",40,0,20,250,0,50); fPhoVertexReco_Invmass->Sumw2(); fOutputList->Add(fPhoVertexReco_Invmass); fPhoVertexReco_step0= new TH2D("fPhoVertexReco_step0","photon production Vertex mass selection",40,0,20,250,0,50); fPhoVertexReco_step0->Sumw2(); fOutputList->Add(fPhoVertexReco_step0); fPhoVertexReco_step1= new TH2D("fPhoVertexReco_step1","photon production Vertex mass selection",40,0,20,250,0,50); fPhoVertexReco_step1->Sumw2(); fOutputList->Add(fPhoVertexReco_step1); fMatchV0_0 = new TH1D("fMatchV0_0","V0 match",100,0,20); fOutputList->Add(fMatchV0_0); fMatchV0_1 = new TH1D("fMatchV0_1","V0 match",100,0,20); fOutputList->Add(fMatchV0_1); fMatchMC_0 = new TH1D("fMatchMC_0","MC match",100,0,20); fOutputList->Add(fMatchMC_0); fMatchMC_1 = new TH1D("fMatchMC_1","MC match",100,0,20); fOutputList->Add(fMatchMC_1); fpair = new TH2D("fpair","pair of associate",100,0,20,21,-10.5,10.5); fOutputList->Add(fpair); fFakeRejection0 = new TH2D("fFakeRejection0","TPC PID",2,-0.5,1.5,100,0,20); fFakeRejection0->Sumw2(); fOutputList->Add(fFakeRejection0); fFakeRejection1 = new TH2D("fFakeRejection1","TPC PID + Tr match",2,-0.5,1.5,100,0,20); fFakeRejection1->Sumw2(); fOutputList->Add(fFakeRejection1); fFakeRejection2 = new TH2D("fFakeRejection2","TPC PID + Tr match + E/p",2,-0.5,1.5,100,0,20); fFakeRejection2->Sumw2(); fOutputList->Add(fFakeRejection2); EopFake = new TH2D("EopFake","negative track Eop",20,0,20,200,0,2); EopFake->Sumw2(); fOutputList->Add(EopFake); EopTrue = new TH2D("EopTrue","true track Eop",20,0,20,200,0,2); EopTrue->Sumw2(); fOutputList->Add(EopTrue); MatchFake = new TH2D("MatchFake","negative track Match",20,0,20,100,0,0.05); MatchFake->Sumw2(); fOutputList->Add(MatchFake); MatchTrue = new TH2D("MatchTrue","true track Match",20,0,20,100,0.05); MatchTrue->Sumw2(); fOutputList->Add(MatchTrue); PostData(1,fOutputList); } //________________________________________________________________________ void AliAnalysisTaskHFECal::Terminate(Option_t *) { // Info("Terminate"); AliAnalysisTaskSE::Terminate(); } //________________________________________________________________________ Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track) { // Check single track cuts for a given cut step const Int_t kMCOffset = AliHFEcuts::kNcutStepsMCTrack; if(!fCFM->CheckParticleCuts(cutStep + kMCOffset, track)) return kFALSE; return kTRUE; } //_________________________________________ void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonic, Bool_t &fFlagConvinat, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent, Bool_t tagpi0, Bool_t tageta, Int_t iCal) { //Identify non-heavy flavour electrons using Invariant mass method fTrackCuts->SetAcceptKinkDaughters(kFALSE); fTrackCuts->SetRequireTPCRefit(kTRUE); fTrackCuts->SetRequireITSRefit(kTRUE); fTrackCuts->SetEtaRange(-0.9,0.9); fTrackCuts->SetMaxChi2PerClusterTPC(3.5); fTrackCuts->SetMinNClustersTPC(90); //const AliESDVertex *pVtx = fESD->GetPrimaryVertex(); Double_t bfield = fESD->GetMagneticField(); Double_t emass = 0.51*0.001; // (0.51 MeV) double ptEle = track->Pt(); //add if(ibgevent==0 && w > 0.0) { fpTCheck->Fill(ptEle,w); } Bool_t flagPhotonicElec = kFALSE; Bool_t flagConvinatElec = kFALSE; int p1 = 0; if(mce==3) { Int_t label = TMath::Abs(track->GetLabel()); TParticle* particle = stack->Particle(label); p1 = particle->GetFirstMother(); } int numULS = 0; int numLS = 0; for(Int_t jTracks = 0; jTracksGetNumberOfTracks(); jTracks++){ AliESDtrack* trackAsso = fESD->GetTrack(jTracks); if (!trackAsso) { printf("ERROR: Could not receive track %d\n", jTracks); continue; } if(itrack==jTracks)continue; int jbgevent = 0; int p2 = 0; if(mce==3) { Int_t label2 = TMath::Abs(trackAsso->GetLabel()); TParticle* particle2 = stack->Particle(label2); Bool_t MChijing_ass = fMC->IsFromBGEvent(label2); if(MChijing_ass)jbgevent =1; if(particle2->GetFirstMother()>-1) p2 = particle2->GetFirstMother(); } Double_t dEdxAsso = -999., ptPrim=-999., ptAsso=-999., openingAngle = -999.; Double_t mass=999., width = -999; Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE; //ptPrim = track->Pt(); ptPrim = ptEle; dEdxAsso = trackAsso->GetTPCsignal(); ptAsso = trackAsso->Pt(); Int_t chargeAsso = trackAsso->Charge(); Int_t charge = track->Charge(); if(ptAsso AcceptTrack(trackAsso)) continue; double fTPCnSigmaAss = fPID->GetPIDResponse() ? fPID->GetPIDResponse()->NumberOfSigmasTPC(trackAsso, AliPID::kElectron) : 1000; //cout << "fTPCnSigmaAss = " << fTPCnSigmaAss << endl; //cout << "fTPCnSigmaAss Cut = " << fMimNsigassCut << endl; if(fTPCnSigmaAss 5) continue; //cout << "fTPCnSigmaAss a.f. cut = " << fTPCnSigmaAss << endl; Int_t fPDGe1 = 11; Int_t fPDGe2 = 11; if(charge>0) fPDGe1 = -11; if(chargeAsso>0) fPDGe2 = -11; //printf("chargeAsso = %d\n",chargeAsso); //printf("charge = %d\n",charge); if(charge == chargeAsso) fFlagLS = kTRUE; if(charge != chargeAsso) fFlagULS = kTRUE; //printf("fFlagLS = %d\n",fFlagLS); //printf("fFlagULS = %d\n",fFlagULS); //printf("\n"); AliKFParticle::SetField(bfield); AliKFParticle ge1(*track, fPDGe1); AliKFParticle ge2(*trackAsso, fPDGe2); if(fSetMassNonlinear) { ge1.SetNonlinearMassConstraint(emass); ge2.SetNonlinearMassConstraint(emass); } AliKFParticle recg(ge1, ge2); /* // vertex AliKFVertex primV(*pVtx); primV += recg; primV -= ge1; primV -= ge2; recg.SetProductionVertex(primV); */ // check chi2 if(recg.GetNDF()<1) continue; Double_t chi2recg = recg.GetChi2()/recg.GetNDF(); Double_t chi2cut = 3.0; // mass. if(fSetMassConstraint) { recg.SetMassConstraint(0,0.0001); chi2cut = 30.0; } recg.GetMass(mass,width); if(fSetMassWidthCut && width>1e10)continue; // angle openingAngle = ge1.GetAngle(ge2); if(fFlagLS) fOpeningAngleLS->Fill(openingAngle); if(fFlagULS) fOpeningAngleULS->Fill(openingAngle); double ishower = 0; if(shower>0.0 && shower<0.3)ishower = 1; if(iCal==1) { double phoinfo[9]; phoinfo[0] = cent; phoinfo[1] = ptPrim; phoinfo[2] = mass; phoinfo[3] = nSig; phoinfo[4] = openingAngle; phoinfo[5] = ishower; phoinfo[6] = ep; phoinfo[7] = mce; phoinfo[8] = ptAsso; if(fFlagLS) fInvmassLS->Fill(phoinfo); if(fFlagULS) fInvmassULS->Fill(phoinfo); if(fFlagLS && ibgevent==0 && jbgevent==0) fInvmassLSmc->Fill(phoinfo,w); if(fFlagULS && ibgevent==0 && jbgevent==0) { fInvmassULSmc->Fill(phoinfo,w); } //printf("fInvmassCut %f\n",fInvmassCut); //printf("openingAngle %f\n",fOpeningAngleCut); } // angle cut if(openingAngle > fOpeningAngleCut) continue; // chi2 cut //if(TMath::Sqrt(TMath::Abs(chi2recg))>chi2cut) continue; if(chi2recg>chi2cut) continue; if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass); if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass); // check double count if(mass0.0 && shower<0.3 && iCal==1)fPhoElecPtM20Mass->Fill(cent,ptPrim); if(mass0.0 && shower<0.3 && iCal==1)fSameElecPtM20Mass->Fill(cent,ptPrim); // flag if(mass0.0 && shower<0.3 && mce>2.1 && iCal==1)fPhoElecPtMCM20Mass->Fill(cent,ptPrim,w); if(mass0.0 && shower<0.3 && mce>2.1 && iCal==1)fSameElecPtMCM20Mass->Fill(cent,ptPrim,w); //printf("MC data\n"); if(w>0.0) { //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl; if(iCal==1) { if(fFlagLS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassLSmc0->Fill(ptPrim,mass); if(fFlagULS && ibgevent==0 && jbgevent==0 && tagpi0) fInvmassULSmc0->Fill(ptPrim,mass); if(fFlagLS && ibgevent==0 && jbgevent==0 && tageta) fInvmassLSmc1->Fill(ptPrim,mass); if(fFlagULS && ibgevent==0 && jbgevent==0 && tageta) fInvmassULSmc1->Fill(ptPrim,mass); if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassLSmc2->Fill(ptPrim,mass); if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tagpi0) fInvmassULSmc2->Fill(ptPrim,mass); if(fFlagLS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassLSmc3->Fill(ptPrim,mass); if(fFlagULS && ibgevent==0 && jbgevent==0 && (p1==p2) && tageta) fInvmassULSmc3->Fill(ptPrim,mass); } } if(mass0 || numLS>0) { int numPair = numULS-numLS; if(iCal==1)fpair->Fill(ptEle,numPair); } fFlagPhotonic = flagPhotonicElec; fFlagConvinat = flagConvinatElec; } //------------------------------------------- void AliAnalysisTaskHFECal::FindMother(TParticle* part, int &label, int &pid) { //int label = 99999; //int pid = 99999; if(part->GetFirstMother()>-1) { label = part->GetFirstMother(); pid = stack->Particle(label)->GetPdgCode(); } //cout << "Find Mother : label = " << label << " ; pid" << pid << endl; } double AliAnalysisTaskHFECal::GetMCweight(double mcPi0pT) { double weight = 1.0; if(mcPi0pT>0.0 && mcPi0pT<5.0) { weight = 0.323*mcPi0pT/(TMath::Exp(-1.6+0.767*mcPi0pT+0.0285*mcPi0pT*mcPi0pT)); } else { weight = 115.0/(0.718*mcPi0pT*TMath::Power(mcPi0pT,3.65)); } return weight; } double AliAnalysisTaskHFECal::GetMCweightEta(double mcEtapT) { double weight = 1.0; weight = 223.3/TMath::Power((TMath::Exp(-0.17*mcEtapT-0.0322*mcEtapT*mcEtapT)+mcEtapT/1.69),5.65); return weight; } //_________________________________________ void AliAnalysisTaskHFECal::FindTriggerClusters() { //cout << "finding trigger patch" << endl; // constants const int nModuleCols = 2; const int nModuleRows = 5; const int nColsFeeModule = 48; const int nRowsFeeModule = 24; const int nColsFaltroModule = 24; const int nRowsFaltroModule = 12; //const int faltroWidthMax = 20; // part 1, trigger extraction ------------------------------------- Int_t globCol, globRow; //Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0, trigInCut=0; Int_t ntimes=0, nTrigChannel=0, nTrigChannelCut=0; //Int_t trigtimes[faltroWidthMax]; Double_t cellTime[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows]; Double_t cellEnergy[nColsFeeModule*nModuleCols][nRowsFeeModule*nModuleRows]; //Double_t fTrigCutLow = 6; //Double_t fTrigCutHigh = 10; Double_t fTimeCutLow = 469e-09; Double_t fTimeCutHigh = 715e-09; AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" ); // erase trigger maps for(Int_t i = 0; i < nColsFaltroModule*nModuleCols; i++ ) { for(Int_t j = 0; j < nRowsFaltroModule*nModuleRows; j++ ) { ftriggersCut[i][j] = 0; ftriggers[i][j] = 0; ftriggersTime[i][j] = 0; } } Int_t iglobCol=0, iglobRow=0; // go through triggers if( fCaloTrigger->GetEntries() > 0 ) { // needs reset fCaloTrigger->Reset(); while( fCaloTrigger->Next() ) { fCaloTrigger->GetPosition( globCol, globRow ); fCaloTrigger->GetNL0Times( ntimes ); /* // no L0s if( ntimes < 1 ) continue; // get precise timings fCaloTrigger->GetL0Times( trigtimes ); trigInCut = 0; for(Int_t i = 0; i < ntimes; i++ ) { // save the first trigger time in channel if( i == 0 || triggersTime[globCol][globRow] > trigtimes[i] ) triggersTime[globCol][globRow] = trigtimes[i]; //printf("trigger times: %d\n",trigtimes[i]); // check if in cut if(trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh ) trigInCut = 1; fTrigTimes->Fill(trigtimes[i]); } */ //L1 analysis from AliAnalysisTaskEMCALTriggerQA Int_t bit = 0; fCaloTrigger->GetTriggerBits(bit); //cout << "bit = " << bit << endl; Int_t ts = 0; fCaloTrigger->GetL1TimeSum(ts); //cout << "ts = " << ts << endl; if (ts > 0)ftriggers[globCol][globRow] = 1; // number of triggered channels in event nTrigChannel++; // ... inside cut if(ts>0 && (bit >> 6 & 0x1)) { iglobCol = globCol; iglobRow = globRow; nTrigChannelCut++; //cout << "ts cut = " << ts << endl; //cout << "globCol = " << globCol << endl; //cout << "globRow = " << globRow << endl; ftriggersCut[globCol][globRow] = 1; } } // calo trigger entries } // has calo trigger entries // part 2 go through the clusters here ----------------------------------- //cout << " part 2 go through the clusters here ----------------------------------- " << endl; Int_t nCluster=0, nCell=0, iCell=0, gCell=0; Short_t cellAddr, nSACell; Int_t mclabel; //Int_t nSACell, iSACell, mclabel; Int_t iSACell; Double_t cellAmp=0, cellTimeT=0, clusterTime=0, efrac=0; Int_t nSupMod, nModule, nIphi, nIeta, iphi, ieta, gphi, geta, feta, fphi; TRefArray *fCaloClusters = new TRefArray(); fESD->GetEMCALClusters( fCaloClusters ); nCluster = fCaloClusters->GetEntries(); // save all cells times if there are clusters if( nCluster > 0 ){ // erase time map for(Int_t i = 0; i < nColsFeeModule*nModuleCols; i++ ){ for(Int_t j = 0; j < nRowsFeeModule*nModuleRows; j++ ){ cellTime[i][j] = 0.; cellEnergy[i][j] = 0.; } } // get cells AliESDCaloCells *fCaloCells = fESD->GetEMCALCells(); //AliVCaloCells fCaloCells = fESD->GetEMCALCells(); nSACell = fCaloCells->GetNumberOfCells(); for(iSACell = 0; iSACell < nSACell; iSACell++ ){ // get the cell info *fCal fCaloCells->GetCell( iSACell, cellAddr, cellAmp, cellTimeT , mclabel, efrac); // get cell position fGeom->GetCellIndex( cellAddr, nSupMod, nModule, nIphi, nIeta ); fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta); // convert co global phi eta gphi = iphi + nRowsFeeModule*(nSupMod/2); geta = ieta + nColsFeeModule*(nSupMod%2); // save cell time and energy cellTime[geta][gphi] = cellTimeT; cellEnergy[geta][gphi] = cellAmp; } } Int_t nClusterTrig, nClusterTrigCut; UShort_t *cellAddrs; Double_t clsE=-999, clsEta=-999, clsPhi=-999; Float_t clsPos[3] = {0.,0.,0.}; for(Int_t icl=0; iclGetNumberOfCaloClusters(); icl++) { AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl); if(!cluster || !cluster->IsEMCAL()) continue; // get cluster cells nCell = cluster->GetNCells(); // get cluster energy clsE = cluster->E(); // get cluster position cluster->GetPosition(clsPos); TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]); clsEta = clsPosVec.Eta(); clsPhi = clsPosVec.Phi(); // get the cell addresses cellAddrs = cluster->GetCellsAbsId(); // check if the cluster contains cell, that was marked as triggered nClusterTrig = 0; nClusterTrigCut = 0; // loop the cells to check, if cluser in acceptance // any cluster with a cell outside acceptance is not considered for( iCell = 0; iCell < nCell; iCell++ ) { // check hot cell //if(clsE>6.0)fCellCheck->Fill(clsE,cellAddrs[iCell]); // get cell position fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta ); fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta); // convert co global phi eta gphi = iphi + nRowsFeeModule*(nSupMod/2); geta = ieta + nColsFeeModule*(nSupMod%2); if( cellTime[geta][gphi] > 0. ){ clusterTime += cellTime[geta][gphi]; gCell++; } // get corresponding FALTRO fphi = gphi / 2; feta = geta / 2; //cout << "fphi = " << fphi << endl; //cout << "feta = " << feta << endl; // try to match with a triggered if( ftriggers[feta][fphi]==1) { nClusterTrig++; } if( ftriggersCut[feta][fphi]==1) { nClusterTrigCut++; } //cout << "nClusterTrigCut : " << nClusterTrigCut << endl; } // cells if( gCell > 0 ) clusterTime = clusterTime / (Double_t)gCell; // fix the reconstriction code time 100ns jumps if( fESD->GetBunchCrossNumber() % 4 < 2 ) clusterTime -= 0.0000001; //fClsETime->Fill(clsE,clusterTime); //fClsEBftTrigCut->Fill(clsE); if(nClusterTrig>0){ //fClsETime1->Fill(clsE,clusterTime); } if(nClusterTrig>0){ cluster->SetChi2(1); //fClsEAftTrigCut1->Fill(clsE); } if(nClusterTrigCut>0){ cluster->SetChi2(2); //fClsEAftTrigCut2->Fill(clsE); } if(nClusterTrigCut>0 && ( nCell > (1 + clsE / 3))) { cluster->SetChi2(3); //fClsEAftTrigCut3->Fill(clsE); } if(nClusterTrigCut>0 && (nCell > (1 + clsE / 3) )&&( clusterTime > fTimeCutLow && clusterTime < fTimeCutHigh )) { // cluster->SetChi2(4); //fClsEAftTrigCut4->Fill(clsE); } if(nClusterTrigCut<1) { cluster->SetChi2(0); //fClsEAftTrigCut->Fill(clsE); } } // clusters } // <-------- only MC correction double AliAnalysisTaskHFECal::MCEopMeanCorrection(double pTmc, float central) { TF1 *fcorr0 = new TF1("fcorr0","[0]*tanh([1]+[2]*x)"); TF1 *fcorr1 = new TF1("fcorr1","[0]*tanh([1]+[2]*x)"); double shift = 0.0; if(central>0 && central<=10) { fcorr0->SetParameters(1.045,1.288,3.18e-01); // fcorr1->SetParameters(9.91e-01,3.466,2.344); } else if(central>10 && central<=20) { fcorr0->SetParameters(1.029,8.254e-01,4.07e-01); fcorr1->SetParameters(0.975,2.276,1.501e-01); } else if(central>20 && central<=30) { fcorr0->SetParameters(1.01,8.795e-01,3.904e-01); fcorr1->SetParameters(9.675e-01,1.654,2.583e-01); } else if(central>30 && central<=40) { fcorr0->SetParameters(1.00,1.466,2.305e-1); fcorr1->SetParameters(9.637e-01,1.473,2.754e-01); } else if(central>40 && central<=50) { fcorr0->SetParameters(1.00,1.422,1.518e-01); fcorr1->SetParameters(9.59e-01,1.421,2.931e-01); } else if(central>50 && central<=70) { fcorr0->SetParameters(0.989,2.495,2.167); fcorr1->SetParameters(0.961,1.734,1.438e-01); } else if(central>70 && central<=100) { fcorr0->SetParameters(0.981,-3.373,3.93327); fcorr1->SetParameters(9.574e-01,1.698,1.58e-01); } shift = fcorr0->Eval(pTmc)-fcorr1->Eval(pTmc); fcorr0->Delete(); fcorr1->Delete(); return shift; } // <-------- only Data correction double AliAnalysisTaskHFECal::NsigmaCorrection(double tmpeta, float central) { int icent = 0; if(central>=0 && central<10) { icent = 0; } else if(central>=10 && central<20) { icent = 1; } else if(central>=20 && central<30) { icent = 2; } else if(central>=30 && central<40) { icent = 3; } else if(central>=40 && central<50) { icent = 4; } else if(central>=50 && central<70) { icent = 5; } else { icent = 6; } double shift = fnSigEtaCorr[icent]->Eval(tmpeta); //cout << "eta correction"<< endl; //cout << "cent = "<< central<< endl; //cout << "icent = "<< icent << endl; //cout << "shift = "<< shift << endl; return shift; }