X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGHF%2Fhfe%2FAliAnalysisTaskHFECal.cxx;h=a69e2de6e0d68c4dac81f36a59358b1c8397e3d4;hb=4ade0749a4a9682748127718c2903dd8a55b1666;hp=60d5de2641293c0ebb60d5a695fc336d13751a2c;hpb=02cd82a00bfa5c1e5843bd3ac3b7e98d81ec2826;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGHF/hfe/AliAnalysisTaskHFECal.cxx b/PWGHF/hfe/AliAnalysisTaskHFECal.cxx index 60d5de26412..a69e2de6e0d 100644 --- a/PWGHF/hfe/AliAnalysisTaskHFECal.cxx +++ b/PWGHF/hfe/AliAnalysisTaskHFECal.cxx @@ -20,12 +20,16 @@ #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" @@ -80,12 +84,15 @@ #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(1) @@ -100,9 +107,14 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) ,fPID(0) ,fPIDqa(0) ,fOpeningAngleCut(0.1) - ,fInvmassCut(0.01) + ,fMimpTassCut(0.5) + ,fInvmassCut(0) // no mass + ,fSetMassConstraint(kTRUE) + ,fSetMassWidthCut(kFALSE) + ,fSetKFpart(kTRUE) ,fNoEvents(0) ,fEMCAccE(0) + ,hEMCAccE(0) ,fTrkpt(0) ,fTrkEovPBef(0) ,fTrkEovPAft(0) @@ -112,6 +124,18 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) ,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) @@ -124,6 +148,7 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) ,fTPCnsigma(0) ,fCent(0) ,fEleInfo(0) + ,fElenSigma(0) /* ,fClsEBftTrigCut(0) ,fClsEAftTrigCut(0) @@ -142,18 +167,54 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal(const char *name) ,fIncpTMChfeAll(0) ,fIncpTMCM20hfe(0) ,fIncpTMCM20hfeAll(0) + ,fIncpTMCM20hfeCheck(0) + ,fInputHFEMC_weight(0) + ,fIncpTMCM20hfeCheck_weight(0) ,fIncpTMCpho(0) ,fIncpTMCM20pho(0) ,fPhoElecPtMC(0) ,fPhoElecPtMCM20(0) ,fSameElecPtMC(0) ,fSameElecPtMCM20(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) + ,CheckNclust(0) + ,CheckNits(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) + //,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()); @@ -169,6 +230,7 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskHFECal") ,fESD(0) ,fMC(0) + ,stack(0) ,fGeom(0) ,fOutputList(0) ,fqahist(1) @@ -183,9 +245,14 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() ,fPID(0) ,fPIDqa(0) ,fOpeningAngleCut(0.1) - ,fInvmassCut(0.01) + ,fMimpTassCut(0.5) + ,fInvmassCut(0) // no mass + ,fSetMassConstraint(kTRUE) + ,fSetMassWidthCut(kFALSE) + ,fSetKFpart(kTRUE) ,fNoEvents(0) ,fEMCAccE(0) + ,hEMCAccE(0) ,fTrkpt(0) ,fTrkEovPBef(0) ,fTrkEovPAft(0) @@ -195,6 +262,18 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() ,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) @@ -207,6 +286,7 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() ,fTPCnsigma(0) ,fCent(0) ,fEleInfo(0) + ,fElenSigma(0) /* ,fClsEBftTrigCut(0) ,fClsEAftTrigCut(0) @@ -225,18 +305,54 @@ AliAnalysisTaskHFECal::AliAnalysisTaskHFECal() ,fIncpTMChfeAll(0) ,fIncpTMCM20hfe(0) ,fIncpTMCM20hfeAll(0) + ,fIncpTMCM20hfeCheck(0) + ,fInputHFEMC_weight(0) + ,fIncpTMCM20hfeCheck_weight(0) ,fIncpTMCpho(0) ,fIncpTMCM20pho(0) ,fPhoElecPtMC(0) ,fPhoElecPtMCM20(0) ,fSameElecPtMC(0) ,fSameElecPtMCM20(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) + ,CheckNclust(0) + ,CheckNits(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) + //,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 @@ -256,6 +372,7 @@ AliAnalysisTaskHFECal::~AliAnalysisTaskHFECal() delete fOutputList; delete fGeom; delete fPID; + delete fCuts; delete fCFM; delete fPIDqa; delete fTrackCuts; @@ -286,7 +403,7 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) } if(fmcData)fMC = MCEvent(); - AliStack* stack = NULL; + //AliStack* stack = NULL; if(fmcData && fMC)stack = fMC->Stack(); Float_t cent = -1.; @@ -308,14 +425,47 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) 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(pTMC,iHijing); + //if(fPDG==221)HETApTcheck->Fill(pTMC,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); + 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); } @@ -332,7 +482,7 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) if(TMath::Abs(pVtxZ)>10) return; fNoEvents->Fill(1); - if(fNOtrks<2) return; + if(fNOtrks<1) return; fNoEvents->Fill(2); AliPIDResponse *pidResponse = fInputHandler->GetPIDResponse(); @@ -357,10 +507,11 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) FindTriggerClusters(); // make EMCAL array + double maxE = 0.0; for(Int_t iCluster=0; iClusterGetNumberOfCaloClusters(); iCluster++) { AliESDCaloCluster *clust = fESD->GetCaloCluster(iCluster); - if(clust->IsEMCAL()) + if(clust && clust->IsEMCAL()) { double clustE = clust->E(); float emcx[3]; // cluster pos @@ -370,8 +521,11 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) double emceta = clustpos.Eta(); double calInfo[5]; calInfo[0] = emcphi; calInfo[1] = emceta; calInfo[2] = clustE; calInfo[3] = cent; calInfo[4] = clust->Chi2(); - //if(clustE>1.5)fEMCAccE->Fill(calInfo); + //fEMCAccE->Fill(calInfo); + //if(clustE>3.0)fEMCAccE->Fill(calInfo); //if(fqahist==1 && clustE>1.5)fEMCAccE->Fill(calInfo); + hEMCAccE->Fill(cent,clustE); + if(clustE>maxE)maxE = clustE; } } @@ -383,32 +537,148 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) continue; } + 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; + + int iHijing = 1; + int mcLabel = -1; + if(fmcData && fMC && stack) { Int_t label = TMath::Abs(track->GetLabel()); - TParticle* particle = stack->Particle(label); - int mcpid = particle->GetPdgCode(); - mcpT = particle->Pt(); - //printf("MCpid = %d",mcpid); - if(particle->GetFirstMother()>-1) - { - int parentPID = stack->Particle(particle->GetFirstMother())->GetPdgCode(); - mcMompT = stack->Particle(particle->GetFirstMother())->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; - } - if(fabs(mcpid)==11 && mcDtoE)mcele= 1.; - if(fabs(mcpid)==11 && mcBtoE)mcele= 2.; - if(fabs(mcpid)==11 && mcPho)mcele= 3.; - } + mcLabel = track->GetLabel(); + + 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(); + //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.0)) continue; @@ -430,6 +700,9 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) // HFE cuts: TPC PID cleanup if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)) continue; + int nTPCcl = track->GetTPCNcls(); + //int nTPCclF = track->GetTPCNclsF(); // warnings + int nITS = track->GetNcls(0); fTrackPtAftTrkCuts->Fill(track->Pt()); @@ -439,14 +712,24 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) // Track extrapolation - Int_t charge = track->Charge(); + //Int_t charge = track->Charge(); fTrkpt->Fill(pt); mom = track->P(); phi = track->Phi(); eta = track->Eta(); 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; + double ncells = -1.0; double m20 = -1.0; double m02 = -1.0; @@ -454,6 +737,8 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) double rmatch = -1.0; double nmatch = -1.0; double oppstatus = 0.0; + double emctof = 0.0; + Bool_t MaxEmatch = kFALSE; Bool_t fFlagPhotonicElec = kFALSE; Bool_t fFlagConvinatElec = kFALSE; @@ -464,7 +749,16 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) if(clust && clust->IsEMCAL()){ double clustE = clust->E(); + if(clustE==maxE)MaxEmatch = kTRUE; eop = clustE/fabs(mom); + //cout << "eop org = "<< eop << endl; + if(mcLabel>-1.0) + { + double mceopcorr = MCEopMeanCorrection(pt,cent); + eop += mceopcorr; + } + //cout << "eop corr = " << eop << endl; + //double clustT = clust->GetTOF(); ncells = clust->GetNCells(); m02 = clust->GetM02(); @@ -474,32 +768,79 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) double deleta = clust->GetTrackDz(); rmatch = sqrt(pow(delphi,2)+pow(deleta,2)); nmatch = clust->GetNTracksMatched(); + emctof = clust->GetTOF(); + //cout << "emctof = " << emctof << endl; if(fTPCnSigma>-1.5 && fTPCnSigma<3.0) { - SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele); + if(fSetKFpart) + { + SelectPhotonicElectron(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta); + } + else + { + SelectPhotonicElectron2(iTracks,cent,track,fFlagPhotonicElec,fFlagConvinatElec,fTPCnSigma,m20,eop,mcele,mcWeight,iHijing,mcOrgPi0,mcOrgEta); + } } if(fFlagPhotonicElec)oppstatus = 1.0; if(fFlagConvinatElec)oppstatus = 2.0; if(fFlagPhotonicElec && fFlagConvinatElec)oppstatus = 3.0; double valdedx[16]; - valdedx[0] = pt; valdedx[1] = dEdx; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma; - valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = m02; valdedx[9] = m20; valdedx[10] = mcpT; - valdedx[11] = cent; valdedx[12] = charge; valdedx[13] = oppstatus; valdedx[14] = clust->Chi2(); + valdedx[0] = pt; valdedx[1] = nITS; valdedx[2] = phi; valdedx[3] = eta; valdedx[4] = fTPCnSigma; + //valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nTPCclF; valdedx[9] = m20; valdedx[10] = mcpT; + valdedx[5] = eop; valdedx[6] = rmatch; valdedx[7] = ncells, valdedx[8] = nmatch; valdedx[9] = m20; valdedx[10] = mcpT; + valdedx[11] = cent; valdedx[12] = dEdx; valdedx[13] = oppstatus; valdedx[14] = nTPCcl; valdedx[15] = mcele; - if(fqahist==1)fEleInfo->Fill(valdedx); + //if(fqahist==1)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); + } + + // ============ PID + + if(nITS<2.5)continue; + if(nTPCcl<100)continue; + + CheckNclust->Fill(nTPCcl); + CheckNits->Fill(nITS); + fdEdxBef->Fill(mom,fTPCnSigma); fTPCnsigma->Fill(mom,fTPCnSigma); if(fTPCnSigma >= -1.0 && fTPCnSigma <= 3)fTrkEovPBef->Fill(pt,eop); Int_t pidpassed = 1; - + + // 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 >= -1.0 && fTPCnSigma <= 3) && mcele>-1 && mcPho && mcOrgPi0) + { + if(iHijing==1)mcWeight = 1.0; + fIncpTMCpho_pi0e_TPC->Fill(phoval,mcWeight); + if(fFlagPhotonicElec) fPhoElecPtMC_pi0e_TPC->Fill(phoval,mcWeight); + if(fFlagConvinatElec) fSameElecPtMC_pi0e_TPC->Fill(phoval,mcWeight); + } //--- track accepted AliHFEpidObject hfetrack; @@ -534,29 +875,70 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) if(m20>0.0 && m20<0.3) { - fIncpTM20->Fill(cent,pt); + fIncpTM20->Fill(cent,pt); + ftimingEle->Fill(pt,emctof); if(fFlagPhotonicElec) fPhoElecPtM20->Fill(cent,pt); if(fFlagConvinatElec) fSameElecPtM20->Fill(cent,pt); } + + //-------- + + double recopT = SumpT(iTracks,track); + + if(m20>0.0 && m20<0.3) + { + if(MaxEmatch)fIncMaxE->Fill(cent,pt); + if(pt>5.0) + { + fIncReco->Fill(cent,recopT); + if(fFlagPhotonicElec) fPhoReco->Fill(cent,recopT); + if(fFlagConvinatElec) fSamReco->Fill(cent,recopT); + if(MaxEmatch) + { + fIncRecoMaxE->Fill(cent,recopT); + if(fFlagPhotonicElec) fPhoRecoMaxE->Fill(cent,recopT); + if(fFlagConvinatElec) fSamRecoMaxE->Fill(cent,recopT); + } + } + } + // MC - if(mcele>0) + // 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) + if(mcBtoE || mcDtoE) // select B->e & D->e { fIncpTMChfe->Fill(cent,pt); - if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt); + //if(m20>0.0 && m20<0.3)fIncpTMCM20hfe->Fill(cent,pt); + //if(m20>0.0 && m20<0.3)fIncpTMCM20hfeCheck->Fill(cent,mcpT); + 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) + + if(mcPho) // select photonic electrons { - double phoval[3]; - phoval[0] = cent; - phoval[1] = pt; - phoval[2] = fTPCnSigma; fIncpTMCpho->Fill(phoval); if(fFlagPhotonicElec) fPhoElecPtMC->Fill(phoval); @@ -567,6 +949,30 @@ void AliAnalysisTaskHFECal::UserExec(Option_t*) 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); + //fIncpTMCM20pho_pi0e->Fill(phoval); // v5-04-02-AN & v5-04-06-AN + //if(fFlagPhotonicElec) fPhoElecPtMCM20_pi0e->Fill(phoval); + //if(fFlagConvinatElec) fSameElecPtMCM20_pi0e->Fill(phoval); + } + if(mcOrgEta) + { + fIncpTMCM20pho_eta->Fill(phoval,mcWeight); + if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval,mcWeight); + if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval,mcWeight); + //fIncpTMCM20pho_eta->Fill(phoval); + //if(fFlagPhotonicElec) fPhoElecPtMCM20_eta->Fill(phoval); + //if(fFlagConvinatElec) fSameElecPtMCM20_eta->Fill(phoval); + } + // --- eta + } } } } @@ -623,7 +1029,23 @@ void AliAnalysisTaskHFECal::UserCreateOutputObjects() 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->SetCutITSpixel(AliHFEextraCuts::kAny); + fCuts->SetCheckITSLayerStatus(kFALSE); + fCuts->SetVertexRange(10.); + fCuts->SetTOFPIDStep(kFALSE); + fCuts->SetPtRange(2, 50); + fCuts->SetMaxImpactParam(3.,3.); + //--------Initialize correction Framework and Cuts fCFM = new AliCFManager; const Int_t kNcutSteps = AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kNcutStepsRecTrack + AliHFEcuts::kNcutStepsDETrack; @@ -650,7 +1072,10 @@ void AliAnalysisTaskHFECal::UserCreateOutputObjects() 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); + //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); @@ -682,16 +1107,62 @@ void AliAnalysisTaskHFECal::UserCreateOutputObjects() fIncpTM20 = new TH2F("fIncpTM20","HFE pid electro vs. centrality with M20",200,0,100,100,0,50); fOutputList->Add(fIncpTM20); - Int_t nBinspho[9] = { 200, 100, 500, 12, 50, 4, 200, 8, 100}; - Double_t minpho[9] = { 0., 0., 0., -2.5, 0, -0.5, 0,-1.5, 0}; - Double_t maxpho[9] = {100., 50., 0.5, 3.5, 1, 3.5, 2, 6.5, 50}; + 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); - fOutputList->Add(fInvmassLS); + 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); - fOutputList->Add(fInvmassULS); + 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); @@ -718,17 +1189,23 @@ void AliAnalysisTaskHFECal::UserCreateOutputObjects() // Make common binning const Double_t kMinP = 0.; - const Double_t kMaxP = 50.; - const Double_t kTPCSigMim = 40.; - const Double_t kTPCSigMax = 140.; + 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) - Int_t nBins[16] = { 250, 200, 60, 20, 100, 300, 50, 40, 200, 200, 250, 200, 3, 5, 10, 8}; - Double_t min[16] = {kMinP, kTPCSigMim, 1.0, -1.0, -6.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, -1.5, -0.5, -0.5, -1.5}; - Double_t max[16] = {kMaxP, kTPCSigMax, 4.0, 1.0, 4.0, 3.0, 0.1, 40, 2.0, 2.0, 50.0, 100, 1.5, 4.5, 9.5, 6.5}; - fEleInfo = new THnSparseD("fEleInfo", "Electron Info; pT [GeV/c]; TPC signal;phi;eta;nSig; E/p;Rmatch;Ncell;M02;M20;mcpT;Centrality;charge;opp;same;trigCond;MCele", 16, nBins, min, max); + Int_t nBins[16] = { 100, 7, 60, 20, 90, 250, 25, 40, 10, 200, 100, 100, 500, 5, 100, 8}; + Double_t min[16] = {kMinP, -0.5, 1.0, -1.0, -5.0, 0, 0, 0, 0.0, 0.0, 0.0, 0, 0, -0.5, 80, -1.5}; + Double_t max[16] = {kMaxP, 6.5, 4.0, 1.0, 4.0, 2.5, 0.05, 40, 10, 2.0, 20.0, 100, 100, 4.5, 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); 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); @@ -781,29 +1258,169 @@ void AliAnalysisTaskHFECal::UserCreateOutputObjects() 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[3] = { 200, 100, 7}; - Double_t minpho2[3] = { 0., 0., -2.5}; - Double_t maxpho2[3] = {100., 50., 4.5}; - - fIncpTMCpho = new THnSparseD("fIncpTMCpho","MC Pho pid electro vs. centrality",3,nBinspho2,minpho2,maxpho2); + 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",3,nBinspho2,minpho2,maxpho2); + 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",3,nBinspho2,minpho2,maxpho2); + 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",3,nBinspho2,minpho2,maxpho2); + fPhoElecPtMCM20 = new THnSparseD("fPhoElecPtMCM20", "MC Pho-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fPhoElecPtMCM20); - fSameElecPtMC = new THnSparseD("fSameElecPtMC", "MC Same-inclusive electron pt",3,nBinspho2,minpho2,maxpho2); + 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",3,nBinspho2,minpho2,maxpho2); + fSameElecPtMCM20 = new THnSparseD("fSameElecPtMCM20", "MC Same-inclusive electron pt with M20",5,nBinspho2,minpho2,maxpho2); fOutputList->Add(fSameElecPtMCM20); + 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); + //------------- + + 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); + /* + 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); PostData(1,fOutputList); } @@ -824,9 +1441,7 @@ Bool_t AliAnalysisTaskHFECal::ProcessCutStep(Int_t cutStep, AliVParticle *track) return kTRUE; } //_________________________________________ -//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec) -//void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig) -void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce) +void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, AliESDtrack *track, Bool_t &fFlagPhotonicElec, Bool_t &fFlagConvinatElec, Double_t nSig, Double_t shower, Double_t ep, Double_t mce, Double_t w, Int_t ibgevent, Bool_t tagpi0, Bool_t tageta) { //Identify non-heavy flavour electrons using Invariant mass method @@ -839,10 +1454,25 @@ void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, fTrackCuts->SetMinNClustersTPC(90); const AliESDVertex *pVtx = fESD->GetPrimaryVertex(); - + Double_t bfield = fESD->GetMagneticField(); + + 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(); + } + //for(Int_t jTracks = itrack+1; jTracksGetNumberOfTracks(); jTracks++){ for(Int_t jTracks = 0; jTracksGetNumberOfTracks(); jTracks++){ AliESDtrack* trackAsso = fESD->GetTrack(jTracks); @@ -851,12 +1481,25 @@ void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, 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 = track->Pt(); + ptPrim = ptEle; dEdxAsso = trackAsso->GetTPCsignal(); ptAsso = trackAsso->Pt(); @@ -864,7 +1507,7 @@ void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, Int_t charge = track->Charge(); - if(ptAsso <0.5) continue; + if(ptAsso AcceptTrack(trackAsso)) continue; if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1 @@ -879,29 +1522,38 @@ void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, //printf("fFlagLS = %d\n",fFlagLS); //printf("fFlagULS = %d\n",fFlagULS); - //printf("\n"); + printf("\n"); + AliKFParticle::SetField(bfield); AliKFParticle ge1(*track, fPDGe1); AliKFParticle ge2(*trackAsso, fPDGe2); AliKFParticle recg(ge1, ge2); - if(recg.GetNDF()<1) continue; - Double_t chi2recg = recg.GetChi2()/recg.GetNDF(); - if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue; - + // vertex AliKFVertex primV(*pVtx); primV += recg; recg.SetProductionVertex(primV); - recg.SetMassConstraint(0,0.0001); - + // 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); - - recg.GetMass(mass,width); - double ishower = 0; if(shower>0.0 && shower<0.3)ishower = 1; @@ -910,6 +1562,7 @@ void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, phoinfo[1] = ptPrim; phoinfo[2] = mass; phoinfo[3] = nSig; + //phoinfo[3] = dEdxAsso; phoinfo[4] = openingAngle; phoinfo[5] = ishower; phoinfo[6] = ep; @@ -918,29 +1571,312 @@ void AliAnalysisTaskHFECal::SelectPhotonicElectron(Int_t itrack, Double_t cent, 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; - - if(masschi2cut) continue; + if(chi2recg>chi2cut) continue; + + if(fFlagLS ) fInvmassLSreco->Fill(ptPrim,mass); + if(fFlagULS) fInvmassULSreco->Fill(ptPrim,mass); + + // for real data + //printf("mce =%f\n",mce); + if(mce<-0.5) // mce==-1. is real + { + //printf("Real data\n"); + if(mass0.0) + { + //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl; + 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(massSetAcceptKinkDaughters(kFALSE); + fTrackCuts->SetRequireTPCRefit(kTRUE); + fTrackCuts->SetRequireITSRefit(kTRUE); + fTrackCuts->SetEtaRange(-0.9,0.9); + //fTrackCuts->SetRequireSigmaToVertex(kTRUE); + fTrackCuts->SetMaxChi2PerClusterTPC(3.5); + fTrackCuts->SetMinNClustersTPC(90); + + //const AliESDVertex *pVtx = fESD->GetPrimaryVertex(); + Double_t eMass = TDatabasePDG::Instance()->GetParticle(11)->Mass(); //Electron mass in GeV + Double_t bfield = fESD->GetMagneticField(); + + 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(); + } + + //for(Int_t jTracks = itrack+1; jTracksGetNumberOfTracks(); jTracks++){ + 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(massGetLabel()); + 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; + Double_t mass=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 <0.5) continue; + if(!fTrackCuts->AcceptTrack(trackAsso)) continue; + if(dEdxAsso <65 || dEdxAsso>100) continue; //11a pass1 + + 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; + + // mass cal 0 + double xt1 = 0.0, xt2 = 0.0; + Double_t p1at[3] = {0,0,0}; + Double_t p2at[3] = {0,0,0}; + //double dca12 = trackAsso->GetDCA(track,bfield,xt2,xt1); //DCA track1-track2 + double kHasdcaT1 = track->GetPxPyPzAt(xt1,bfield,p1at); //Track1 + double kHasdcaT2 = trackAsso->GetPxPyPzAt(xt2,bfield,p2at); //Track2 + if(!kHasdcaT1 || !kHasdcaT2) AliWarning("It could be a problem in the extrapolation"); + //cout << "dca = " << dca12 << endl; + // 3D + TLorentzVector electron1; + TLorentzVector electron2; + TLorentzVector mother; + + electron1.SetXYZM(p1at[0], p1at[1], p1at[2], eMass); + electron2.SetXYZM(p2at[0], p2at[1], p2at[2], eMass); + openingAngle = TVector2::Phi_0_2pi(electron1.Angle(electron2.Vect())); + + mother = electron1 + electron2; + //double invmassAtDCA = mother.M(); + + // 2D + TLorentzVector electron1_2D; + TLorentzVector electron2_2D; + TLorentzVector mother_2D; + + double pT1at = sqrt(pow(p1at[0],2)+pow(p1at[1],2)); + double pT2at = sqrt(pow(p2at[0],2)+pow(p2at[1],2)); + + electron1_2D.SetXYZM(pT1at, 0.0, p1at[2], eMass); + electron2_2D.SetXYZM(pT2at, 0.0, p2at[2], eMass); + + mother_2D = electron1_2D + electron2_2D; + double invmassAtDCA_2D = mother_2D.M(); + mass = invmassAtDCA_2D; + + // angle + if(fFlagLS) fOpeningAngleLS->Fill(openingAngle); + if(fFlagULS) fOpeningAngleULS->Fill(openingAngle); + + double ishower = 0; + if(shower>0.0 && shower<0.3)ishower = 1; + + double phoinfo[9]; + phoinfo[0] = cent; + phoinfo[1] = ptPrim; + phoinfo[2] = mass; + phoinfo[3] = nSig; + //phoinfo[3] = dEdxAsso; + 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); + + // for real data + //printf("mce =%f\n",mce); + if(mce<-0.5) // mce==-1. is real + { + //printf("Real data\n"); + if(mass0.0) + { + //cout << "tagpi0 = " << tagpi0 << " ; tageta = " << tageta << endl; + 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(massGetFirstMother()>-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; @@ -1009,9 +1945,11 @@ void AliAnalysisTaskHFECal::FindTriggerClusters() //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++; @@ -1021,6 +1959,9 @@ void AliAnalysisTaskHFECal::FindTriggerClusters() iglobCol = globCol; iglobRow = globRow; nTrigChannelCut++; + //cout << "ts cut = " << ts << endl; + //cout << "globCol = " << globCol << endl; + //cout << "globRow = " << globRow << endl; ftriggersCut[globCol][globRow] = 1; } @@ -1028,8 +1969,10 @@ void AliAnalysisTaskHFECal::FindTriggerClusters() } // 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, mclabel; + 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; @@ -1126,6 +2069,9 @@ void AliAnalysisTaskHFECal::FindTriggerClusters() 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++; @@ -1134,6 +2080,8 @@ void AliAnalysisTaskHFECal::FindTriggerClusters() { nClusterTrigCut++; } + //cout << "nClusterTrigCut : " << nClusterTrigCut << endl; + } // cells @@ -1181,9 +2129,138 @@ void AliAnalysisTaskHFECal::FindTriggerClusters() } // 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); + 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; + +} + + +double AliAnalysisTaskHFECal::SumpT(Int_t itrack, AliESDtrack* track) +{ + + fTrackCuts->SetAcceptKinkDaughters(kFALSE); + fTrackCuts->SetRequireTPCRefit(kTRUE); + fTrackCuts->SetRequireITSRefit(kTRUE); + fTrackCuts->SetEtaRange(-0.9,0.9); + //fTrackCuts->SetRequireSigmaToVertex(kTRUE); + fTrackCuts->SetMaxChi2PerClusterTPC(3.5); + fTrackCuts->SetMinNClustersTPC(90); + + double pTrecp = track->Pt(); + double phiorg = track->Phi(); + double etaorg = track->Eta(); + + 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; + double pTAss = trackAsso->Pt(); + double etaAss = trackAsso->Eta(); + double phiAss = trackAsso->Phi(); + + double delphi = phiorg - phiAss; + double deleta = etaorg - etaAss; + + double R = sqrt(pow(deleta,2)+pow(delphi,2)); + if(pTAss<0.5)continue; + if(R<0.4)pTrecp+=pTAss; + + } + + return pTrecp; +}