X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGGA%2FEMCALTasks%2FAliAnalysisTaskEMCALIsoPhoton.cxx;h=9ac6a1c17b2406f6432ade269ac5831711985843;hb=abed0cc7329c0b7819fcd92fc87e2c1af933468f;hp=d6c46fa5bf955dc25b41f7f45debc89861d61758;hpb=672df96f0184532d3e9f74d1350e7453bc57877c;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx index d6c46fa5bf9..9ac6a1c17b2 100644 --- a/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx +++ b/PWGGA/EMCALTasks/AliAnalysisTaskEMCALIsoPhoton.cxx @@ -24,24 +24,36 @@ #include "AliESDUtils.h" #include "AliESDtrack.h" #include "AliESDtrackCuts.h" +#include "AliAODEvent.h" +#include "AliAODTrack.h" #include "AliMCEvent.h" #include "AliMCEventHandler.h" #include "AliStack.h" +#include "AliVEvent.h" +#include "AliVTrack.h" #include "AliV0vertexer.h" #include "AliVCluster.h" +#include "AliOADBContainer.h" + +#include +using std::cout; +using std::endl; ClassImp(AliAnalysisTaskEMCALIsoPhoton) //________________________________________________________________________ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : AliAnalysisTaskSE(), - fCaloClusters(0), + fESDClusters(0), + fAODClusters(0), fSelPrimTracks(0), fTracks(0), - fEMCalCells(0), + fESDCells(0), + fAODCells(0), fPrTrCuts(0), fGeom(0x0), fGeoName("EMCAL_COMPLETEV1"), + fOADBContainer(0), fPeriod("LHC11c"), fTrigBit("kEMC7"), fIsTrain(0), @@ -57,7 +69,13 @@ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : fNClusForDirPho(0), fDirPhoPt(0), fHigherPtCone(0), + fImportGeometryFromFile(0), + fImportGeometryFilePath(""), + fMaxPtTrack(0), + fMaxEClus(0), + fNCells50(0), fESD(0), + fAOD(0), fMCEvent(0), fStack(0), fOutputList(0), @@ -67,6 +85,7 @@ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : fPVtxZ(0), fTrMultDist(0), fMCDirPhotonPtEtaPhi(0), + fMCIsoDirPhotonPtEtaPhi(0), fDecayPhotonPtMC(0), fCellAbsIdVsAmpl(0), fNClusHighClusE(0), @@ -81,21 +100,41 @@ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : fAllIsoEtMcGamma(0), fAllIsoNoUeEtMcGamma(0), fMCDirPhotonPtEtaPhiNoClus(0), - fHnOutput(0) + fHnOutput(0), + fQAList(0), + fNTracks(0), + fEmcNCells(0), + fEmcNClus(0), + fEmcNClusCut(0), + fNTracksECut(0), + fEmcNCellsCut(0), + fEmcClusEPhi(0), + fEmcClusEPhiCut(0), + fEmcClusEEta(0), + fEmcClusEEtaCut(0), + fTrackPtPhi(0), + fTrackPtPhiCut(0), + fTrackPtEta(0), + fTrackPtEtaCut(0), + fMaxCellEPhi(0) { // Default constructor. + for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; } //________________________________________________________________________ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : AliAnalysisTaskSE(name), - fCaloClusters(0), + fESDClusters(0), + fAODClusters(0), fSelPrimTracks(0), fTracks(0), - fEMCalCells(0), + fESDCells(0), + fAODCells(0), fPrTrCuts(0), fGeom(0x0), fGeoName("EMCAL_COMPLETEV1"), + fOADBContainer(0), fPeriod("LHC11c"), fTrigBit("kEMC7"), fIsTrain(0), @@ -111,7 +150,13 @@ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : fNClusForDirPho(0), fDirPhoPt(0), fHigherPtCone(0), + fImportGeometryFromFile(0), + fImportGeometryFilePath(""), + fMaxPtTrack(0), + fMaxEClus(0), + fNCells50(0), fESD(0), + fAOD(0), fMCEvent(0), fStack(0), fOutputList(0), @@ -121,6 +166,7 @@ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : fPVtxZ(0), fTrMultDist(0), fMCDirPhotonPtEtaPhi(0), + fMCIsoDirPhotonPtEtaPhi(0), fDecayPhotonPtMC(0), fCellAbsIdVsAmpl(0), fNClusHighClusE(0), @@ -135,7 +181,23 @@ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : fAllIsoEtMcGamma(0), fAllIsoNoUeEtMcGamma(0), fMCDirPhotonPtEtaPhiNoClus(0), - fHnOutput(0) + fHnOutput(0), + fQAList(0), + fNTracks(0), + fEmcNCells(0), + fEmcNClus(0), + fEmcNClusCut(0), + fNTracksECut(0), + fEmcNCellsCut(0), + fEmcClusEPhi(0), + fEmcClusEPhiCut(0), + fEmcClusEEta(0), + fEmcClusEEtaCut(0), + fTrackPtPhi(0), + fTrackPtPhiCut(0), + fTrackPtEta(0), + fTrackPtEtaCut(0), + fMaxCellEPhi(0) { // Constructor @@ -145,6 +207,7 @@ AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : // Output slot #0 id reserved by the base class for AOD // Output slot #1 writes into a TH1 container DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); } //________________________________________________________________________ @@ -152,7 +215,7 @@ void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects() { // Create histograms, called once. - fCaloClusters = new TRefArray(); + fESDClusters = new TObjArray(); fSelPrimTracks = new TObjArray(); @@ -160,7 +223,9 @@ void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects() fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data()); - + fOADBContainer = new AliOADBContainer("AliEMCALgeo"); + fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo"); + fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2); fOutputList->Add(fEvtSel); @@ -176,10 +241,14 @@ void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects() fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5); fOutputList->Add(fTrMultDist); - fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20); + fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",100,-0.5,99.5,154,-0.77,0.77,130,1.38,3.20); fMCDirPhotonPtEtaPhi->Sumw2(); fOutputList->Add(fMCDirPhotonPtEtaPhi); + fMCIsoDirPhotonPtEtaPhi = new TH3F("hMCIsoDirPhotonPtEtaPhi","photon (gq->#gammaq, isolated@MC) p_{T}, #eta, #phi;GeV/c;#eta;#phi",100,-0.5,99.5,154,-0.77,0.77,130,1.38,3.20); + fMCIsoDirPhotonPtEtaPhi->Sumw2(); + fOutputList->Add(fMCIsoDirPhotonPtEtaPhi); + fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",1000,0,100); fDecayPhotonPtMC->Sumw2(); fOutputList->Add(fDecayPhotonPtMC); @@ -224,25 +293,83 @@ void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects() fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",1000,0,100,154,-0.77,0.77,130,1.38,3.20); fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus); - Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000, nAllIso=1000, nCeIsoNoUE=1000, nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=100, nPhoMcPt=101; + Int_t nEt=1000, nM02=400, nCeIso=1000, nTrIso=1000, nAllIso=1000, nCeIsoNoUE=1000, nAllIsoNoUE=1000, nTrClDphi=200, nTrClDeta=100, nClEta=140, nClPhi=128, nTime=60, nMult=100, nPhoMcPt=100; Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt}; fNDimensions = sizeof(bins)/sizeof(Int_t); const Int_t ndims = fNDimensions; - Double_t xmin[] = { 0., 0., -10., -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-1.5}; - Double_t xmax[] = { 100., 4., 190., 190., 190., 190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5}; + Double_t xmin[] = { -0.5, 0., -10., -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,-0.5}; + Double_t xmax[] = { 99.5, 4., 190., 190., 190., 190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,99.5}; + if(fPeriod.Contains("11h")){ + xmax[12]=3999.5; + } fHnOutput = new THnSparseF("fHnOutput","Output matrix: E_{T},M02,CeIso,TrIso,AllIso, CeIsoNoUESub, AllIsoNoUESub, d#phi_{trk},d#eta_{trk},#eta_{clus},#phi_{clus},T_{max},mult,mc-p_{T}^{#gamma}", ndims, bins, xmin, xmax); fOutputList->Add(fHnOutput); + //QA outputs + fQAList = new TList(); + fQAList->SetOwner();// Container cleans up all histos (avoids leaks in merging) + + fNTracks = new TH1F("hNTracks","# of selected tracks;n-tracks;counts",120,-0.5,119.5); + fNTracks->Sumw2(); + fQAList->Add(fNTracks); + + fEmcNCells = new TH1F("fEmcNCells",";n/event;count",120,-0.5,119.5); + fEmcNCells->Sumw2(); + fQAList->Add(fEmcNCells); + fEmcNClus = new TH1F("fEmcNClus",";n/event;count",120,-0.5,119.5); + fEmcNClus->Sumw2(); + fQAList->Add(fEmcNClus); + fEmcNClusCut = new TH1F("fEmcNClusCut",";n/event;count",120,-0.5,119.5); + fEmcNClusCut->Sumw2(); + fQAList->Add(fEmcNClusCut); + fNTracksECut = new TH1F("fNTracksECut",";n/event;count",120,-0.5,119.5); + fNTracksECut->Sumw2(); + fQAList->Add(fNTracksECut); + fEmcNCellsCut = new TH1F("fEmcNCellsCut",";n/event;count",120,-0.5,119.5); + fEmcNCellsCut->Sumw2(); + fQAList->Add(fEmcNCellsCut); + fEmcClusEPhi = new TH2F("fEmcClusEPhi",";GeV;#phi",100,-0.25,49.75,63,0,6.3); + fEmcClusEPhi->Sumw2(); + fQAList->Add(fEmcClusEPhi); + fEmcClusEPhiCut = new TH2F("fEmcClusEPhiCut",";GeV;#phi",100,-0.25,49.75,63,0,6.3); + fEmcClusEPhiCut->Sumw2(); + fQAList->Add(fEmcClusEPhiCut); + fEmcClusEEta = new TH2F("fEmcClusEEta",";GeV;#eta",100,-0.25,49.75,19,-0.9,0.9); + fEmcClusEEta->Sumw2(); + fQAList->Add(fEmcClusEEta); + fEmcClusEEtaCut = new TH2F("fEmcClusEEtaCut",";GeV;#eta",100,-0.25,49.75,18,-0.9,0.9); + fEmcClusEEtaCut->Sumw2(); + fQAList->Add(fEmcClusEEtaCut); + fTrackPtPhi = new TH2F("fTrackPtPhi",";p_{T} [GeV/c];#phi",100,-0.25,49.75,63,0,6.3); + fTrackPtPhi->Sumw2(); + fQAList->Add(fTrackPtPhi); + fTrackPtPhiCut = new TH2F("fTrackPtPhiCut",";p_{T} [GeV/c];#phi",100,-0.25,49.75,63,0,6.3); + fTrackPtPhiCut->Sumw2(); + fQAList->Add(fTrackPtPhiCut); + fTrackPtEta = new TH2F("fTrackPtEta",";p_{T} [GeV/c];#eta",100,-0.25,49.75,18,-0.9,0.9); + fTrackPtEta->Sumw2(); + fQAList->Add(fTrackPtEta); + fTrackPtEtaCut = new TH2F("fTrackPtEtaCut",";p_{T} [GeV/c];#eta",100,-0.25,49.75,18,-0.9,0.9); + fTrackPtEtaCut->Sumw2(); + fQAList->Add(fTrackPtEtaCut); + + fMaxCellEPhi = new TH2F("fMaxCellEPhi","Most energetic cell in event; GeV;#phi",100,-0.25,49.75,63,0,6.3); + fMaxCellEPhi->Sumw2(); + fQAList->Add(fMaxCellEPhi); PostData(1, fOutputList); + PostData(2, fQAList); } //________________________________________________________________________ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) { // Main loop, called for each event. - + fESDClusters = 0; + fESDCells = 0; + fAODClusters = 0; + fAODCells = 0; // event trigger selection Bool_t isSelected = 0; if(fPeriod.Contains("11a")){ @@ -260,6 +387,12 @@ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) else isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7); } + if(fPeriod.Contains("11h")){ + if(fTrigBit.Contains("kAny")) + isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny); + if(fTrigBit.Contains("kAnyINT")) + isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT); + } if(fIsMc) isSelected = kTRUE; if(fDebug) @@ -273,20 +406,40 @@ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) if(!filename.Contains(fPathStrOpt.Data())) return; } - fESD = dynamic_cast(InputEvent()); - if (!fESD) { - printf("ERROR: fESD not available\n"); + AliVEvent *event = (AliVEvent*)InputEvent(); + if (!event) { + printf("ERROR: event not available\n"); return; } + Int_t runnumber = InputEvent()->GetRunNumber() ; + if(fDebug) + printf("run number = %d\n",runnumber); + + fESD = dynamic_cast(event); + if(!fESD){ + fAOD = dynamic_cast(event); + if(!fAOD){ + printf("ERROR: Invalid type of event!!!\n"); + return; + } + else if(fDebug) + printf("AOD EVENT!\n"); + } fEvtSel->Fill(0); if(fDebug) - printf("fESD is ok\n"); + printf("event is ok,\n run number=%d\n",runnumber); + - AliESDVertex *pv = (AliESDVertex*)fESD->GetPrimaryVertex(); + AliVVertex *pv = (AliVVertex*)event->GetPrimaryVertex(); + Bool_t pvStatus = kTRUE; + if(fESD){ + AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex(); + pvStatus = esdv->GetStatus(); + } if(!pv) return; - if(!pv->GetStatus()) + if(!pvStatus) fRecoPV->Fill(0); else fRecoPV->Fill(1); @@ -297,10 +450,15 @@ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) printf("passed vertex cut\n"); fEvtSel->Fill(1); + if(fESD) + fTracks = dynamic_cast(InputEvent()->FindListObject("Tracks")); + if(fAOD) + fTracks = dynamic_cast(fAOD->GetTracks()); - fTracks = dynamic_cast(InputEvent()->FindListObject("Tracks")); if(!fTracks){ AliError("Track array in event is NULL!"); + if(fDebug) + printf("returning due to not finding tracks!\n"); return; } // Track loop to fill a pT spectrum @@ -308,29 +466,62 @@ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { // for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks); - AliVTrack *track = static_cast(fTracks->At(iTracks)); + AliVTrack *track = (AliVTrack*)fTracks->At(iTracks); if (!track) continue; - if (fPrTrCuts && fPrTrCuts->IsSelected(track)){ + AliAODTrack *aodTrack = dynamic_cast(track); + AliESDtrack *esdTrack = dynamic_cast(track); + if (esdTrack && fPrTrCuts && fPrTrCuts->IsSelected(track)){ fSelPrimTracks->Add(track); + /*if(fTrackMaxPtPt()) + fTrackMaxPt = track->Pt();*/ //printf("pt,eta,phi:%1.1f,%1.1f,%1.1f \n",track->Pt(),track->Eta(), track->Phi()); } + else if(aodTrack){ + fSelPrimTracks->Add(track); + /*if(fTrackMaxPtPt()) + fTrackMaxPt = track->Pt();*/ + } } - if(!fIsTrain){ - for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ - if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3) - break; + TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices"); + for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ + if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3) + break; + /*if(fESD) fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod); - } + else*/ + // if(event->GetEMCALMatrix(mod)) + fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod); + fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod); + } + + if(fESD){ + AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts(); + fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); + if(fDebug) + printf("ESD Track mult= %d\n",fTrackMult); + } + else if(fAOD){ + fTrackMult = Ntracks; + if(fDebug) + printf("AOD Track mult= %d\n",fTrackMult); } - AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts(); - fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); fTrMultDist->Fill(fTrackMult); - fESD->GetEMCALClusters(fCaloClusters); - fEMCalCells = fESD->GetEMCALCells(); - + if(fESD){ + TList *l = fESD->GetList(); + fESDClusters = dynamic_cast(l->FindObject("CaloClusters")); + fESDCells = fESD->GetEMCALCells(); + if(fDebug) + printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast()); + } + else if(fAOD){ + fAODClusters = dynamic_cast(fAOD->GetCaloClusters()); + fAODCells = fAOD->GetEMCALCells(); + if(fDebug) + printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast()); + } fMCEvent = MCEvent(); @@ -343,6 +534,7 @@ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) if(fDebug) std::cout<<"ERROR: NO MC EVENT!!!!!!\n"; } + LoopOnCells(); FollowGamma(); if(fDebug) printf("passed calling of FollowGamma\n"); @@ -352,22 +544,41 @@ void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) FillMcHists(); if(fDebug) printf("passed calling of FillMcHists\n"); - - fCaloClusters->Clear(); + FillQA(); + if(fDebug) + printf("passed calling of FillQA\n"); + /*if(fESD) + fESDClusters->Clear();*/ fSelPrimTracks->Clear(); fNClusForDirPho = 0; + fNCells50 = 0; PostData(1, fOutputList); + PostData(2, fQAList); } //________________________________________________________________________ void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() { + if(fDebug) + printf("Inside FillClusHists()....\n"); // Fill cluster histograms. + TObjArray *clusters = fESDClusters; - if(!fCaloClusters) + if (!clusters){ + clusters = fAODClusters; + if(fDebug) + printf("ESD clusters empty..."); + } + if (!clusters){ + if(fDebug) + printf(" and AOD clusters as well!!!\n"); return; - const Int_t nclus = fCaloClusters->GetEntries(); + } + if(fDebug) + printf("\n"); + + const Int_t nclus = clusters->GetEntries(); if(nclus==0) return; if(fDebug) @@ -377,7 +588,7 @@ void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() Double_t ptmc=-1; for(Int_t ic=0;ic(fCaloClusters->At(ic)); + AliVCluster *c = static_cast(clusters->At(ic)); if(!c) continue; if(!c->IsEMCAL()) @@ -423,8 +634,7 @@ void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() fMcPtInConeBG->Fill(alliso-Et-allisoue, mcptsum); fMcPtInConeBGnoUE->Fill(alliso-Et, mcptsum); } - Double_t r = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx()+c->GetTrackDz()*c->GetTrackDz()); - if(c->GetM02()>0.1 && c->GetM02()<0.3 && r>0.03){ + if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){ fMcPtInConeSBG->Fill(alliso-Et-allisoue, mcptsum); fMcPtInConeSBGnoUE->Fill(alliso-Et, mcptsum); if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){ @@ -434,19 +644,25 @@ void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() } const Int_t ndims = fNDimensions; Double_t outputValues[ndims]; - ptmc = GetClusSource(c); + if(mcptsum<2) + ptmc = GetClusSource(c); + else + ptmc = -0.1; outputValues[0] = Et; outputValues[1] = c->GetM02(); - outputValues[2] = ceiso-Et/*cecore*/-ceisoue; + outputValues[2] = ceiso/*cecore*/-ceisoue; outputValues[3] = triso-trisoue; - outputValues[4] = alliso-Et/*cecore*/-allisoue; - outputValues[5] = ceiso-Et; - outputValues[6] = alliso-Et; + outputValues[4] = alliso/*cecore*/-allisoue - trcore; + outputValues[5] = ceiso; + outputValues[6] = alliso - trcore; outputValues[7] = c->GetTrackDx(); outputValues[8] = c->GetTrackDz(); outputValues[9] = clsVec.Eta(); outputValues[10] = clsVec.Phi(); - outputValues[11] = fEMCalCells->GetCellTime(id); + if(fESDCells) + outputValues[11] = fESDCells->GetCellTime(id); + else if(fAODCells) + outputValues[11] = fAODCells->GetCellTime(id); outputValues[12] = fTrackMult; outputValues[13] = ptmc; fHnOutput->Fill(outputValues); @@ -454,6 +670,7 @@ void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() maxE = c->E(); } fNClusHighClusE->Fill(maxE,nclus); + fMaxEClus = maxE; fNClusEt10->Fill(nclus10); fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho); } @@ -462,24 +679,40 @@ void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core) { // Get cell isolation. + AliVCaloCells *cells = fESDCells; + if (!cells){ + cells = fAODCells; + if(fDebug) + printf("ESD cells empty..."); + } + if (!cells){ + if(fDebug) + printf(" and AOD cells are empty as well!!!\n"); + return; + } - if(!fEMCalCells) + TObjArray *clusters = fESDClusters; + if (!clusters) + clusters = fAODClusters; + if (!clusters) return; - const Int_t ncells = fEMCalCells->GetNumberOfCells(); + + + const Int_t nclus = clusters->GetEntries(); + //const Int_t ncells = cells->GetNumberOfCells(); Float_t totiso=0; Float_t totband=0; Float_t totcore=0; Float_t etacl = vec.Eta(); Float_t phicl = vec.Phi(); - Float_t thetacl = vec.Theta(); - Float_t maxtcl = fEMCalCells->GetCellTime(maxid); + Float_t maxtcl = cells->GetCellTime(maxid); if(phicl<0) phicl+=TMath::TwoPi(); - Int_t absid = -1; + /*Int_t absid = -1; Float_t eta=-1, phi=-1; for(int icell=0;icellGetCellNumber(icell)); - Float_t celltime = fEMCalCells->GetCellTime(absid); + absid = TMath::Abs(cells->GetCellNumber(icell)); + Float_t celltime = cells->GetCellTime(absid); //if(TMath::Abs(celltime)>2e-8 && (!fIsMc)) if(TMath::Abs(celltime-maxtcl)>2e-8 ) continue; @@ -488,26 +721,50 @@ void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t fGeom->EtaPhiFromIndex(absid,eta,phi); Float_t dphi = TMath::Abs(phi-phicl); Float_t deta = TMath::Abs(eta-etacl); + Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/ + for(int ic=0;ic(clusters->At(ic)); + if(!c) + continue; + if(!c->IsEMCAL()) + continue; + Short_t id; + GetMaxCellEnergy( c, id); + Double_t maxct = cells->GetCellTime(id); + if(TMath::Abs(maxtcl-maxct)>2.5e-9) + continue; + Float_t clsPos[3] = {0,0,0}; + c->GetPosition(clsPos); + TVector3 cv(clsPos); + Double_t Et = c->E()*TMath::Sin(cv.Theta()); + Float_t dphi = TMath::Abs(cv.Phi()-phicl); + Float_t deta = TMath::Abs(cv.Eta()-etacl); Float_t R = TMath::Sqrt(deta*deta + dphi*dphi); - Float_t etcell = fEMCalCells->GetCellAmplitude(absid)*TMath::Sin(thetacl); + if(R<0.007) + continue; + if(maxid==id) + continue; + Double_t matchedpt = GetTrackMatchedPt(c->GetTrackMatchedIndex()); + Double_t nEt = TMath::Max(Et-matchedpt, 0.0); + if(nEt<0) + printf("nEt=%1.1f\n",nEt); if(RfIsoConeR) continue; if(detaGetCellAmplitude(absid)*TMath::Sin(thetacl); + totband += nEt; } } iso = totiso; phiband = totband; core = totcore; } - //________________________________________________________________________ void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core) { @@ -558,7 +815,9 @@ Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluste // Calculate the energy of cross cells around the leading cell. AliVCaloCells *cells = 0; - cells = fESD->GetEMCALCells(); + cells = fESDCells; + if (!cells) + cells = fAODCells; if (!cells) return 0; @@ -607,8 +866,10 @@ Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *clu id = -1; AliVCaloCells *cells = 0; - cells = fESD->GetEMCALCells(); + cells = fESDCells; if (!cells) + cells = fAODCells; + if(!cells) return 0; Double_t maxe = 0; @@ -649,6 +910,9 @@ void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists() Int_t pdgMom = mcmom->GetPdgCode(); if((imom==6 || imom==7) && pdgMom==22) { fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); + Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR); + if(ptsum<2) + fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); if(fNClusForDirPho==0) fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); if(fDebug){ @@ -666,15 +930,15 @@ void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists() Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c) { if(!c) - return -1; + return -0.1; if(!fStack) - return -1; + return -0.1; Int_t nmcp = fStack->GetNtrack(); Int_t clabel = c->GetLabel(); if(fDebug && fMcIdFamily.Contains(Form("%d",clabel))) printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel); if(!fMcIdFamily.Contains(Form("%d",clabel))) - return -1; + return -0.1; fNClusForDirPho++; TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1); Int_t partonpos = partonposstr.Atoi(); @@ -688,16 +952,16 @@ Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c) Double_t Et = c->E()*TMath::Sin(clsVec.Theta()); TParticle *mcp = static_cast(fStack->Particle(partonpos)); if(!mcp) - return -1; + return -0.1; if(fDebug){ printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),clabel); } Int_t lab1 = mcp->GetFirstDaughter(); if(lab1<0 || lab1>nmcp) - return -1; + return -0.1; TParticle *mcd = static_cast(fStack->Particle(lab1)); if(!mcd) - return -1; + return -0.1; if(fDebug) printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->Energy(),mcd->GetPdgCode(),lab1); if(mcd->GetPdgCode()==22){ @@ -706,7 +970,7 @@ Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c) } else{ printf("Warning: daughter of photon parton is not a photon\n"); - return -1; + return -0.1; } return fDirPhoPt; } @@ -781,18 +1045,15 @@ Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t Int_t nTracks = fStack->GetNtrack(); Float_t ptsum = 0; TString addpartlabels = ""; - for(Int_t iTrack=0;iTrack(fStack->Particle(iTrack)); if(!mcp) continue; - Int_t firstd = mcp->GetFirstDaughter(); - Int_t lastd = mcp->GetLastDaughter(); - if((iTrack==6 || iTrack==7) && mcp->GetPdgCode()==22){ - for(Int_t id=firstd;id<=lastd;id++) - addpartlabels += Form("%d,",id); - continue; - } - if(mcp->Rho()>2.5) + Int_t status = mcp->GetStatusCode(); + if(status!=1) + continue; + Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz()); + if(mcvr>1) continue; else { if(fDebug) @@ -802,24 +1063,122 @@ Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t Float_t deta = mcp->Eta() - etaclus; if(fDebug) printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta); - if(deta>10) + if(deta>R || dphi>R) continue; Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta); if(dR>R) continue; - if(addpartlabels.Contains(Form("%d",iTrack))){ - if(fDebug) - printf(" >>>> list of particles included (and their daughters) %s\n",addpartlabels.Data()); - continue; - } addpartlabels += Form("%d,",iTrack); - for(Int_t id=firstd;id<=lastd;id++) - addpartlabels += Form("%d,",id); + if(mcp->Pt()<0.2) + continue; ptsum += mcp->Pt(); } return ptsum; } //________________________________________________________________________ +void AliAnalysisTaskEMCALIsoPhoton::FillQA() +{ + if(!fSelPrimTracks) + return; + const int ntracks = fSelPrimTracks->GetEntriesFast(); + const int ncells = fNCells50;//fESDCells->GetNumberOfCells(); + const Int_t nclus = fESDClusters->GetEntries(); + + fNTracks->Fill(ntracks); + fEmcNCells->Fill(ncells); + fEmcNClus->Fill(nclus); + if(fMaxEClus>fECut){ + fNTracksECut->Fill(ntracks); + fEmcNCellsCut->Fill(ncells); + fEmcNClusCut->Fill(nclus); + } + for(int it=0;itAt(it); + if(!t) + continue; + fTrackPtPhi->Fill(t->Pt(),t->Phi()); + fTrackPtEta->Fill(t->Pt(),t->Eta()); + if(fMaxEClus>fECut){ + fTrackPtPhiCut->Fill(t->Pt(), t->Phi()); + fTrackPtEtaCut->Fill(t->Pt(), t->Eta()); + } + } + for(int ic=0;ic(fESDClusters->At(ic)); + //AliESDCaloCluster *c = (AliESDCaloCluster*)fESDClusters->At(ic); + if(!c) + continue; + if(!c->IsEMCAL()) + continue; + Float_t clsPos[3] = {0,0,0}; + c->GetPosition(clsPos); + TVector3 clsVec(clsPos); + Double_t cphi = clsVec.Phi(); + Double_t ceta = clsVec.Eta(); + fEmcClusEPhi->Fill(c->E(), cphi); + fEmcClusEEta->Fill(c->E(), ceta); + if(fMaxEClus>fECut){ + fEmcClusEPhiCut->Fill(c->E(), cphi); + fEmcClusEEtaCut->Fill(c->E(), ceta); + } + } +} +//________________________________________________________________________ +Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex) +{ + Double_t pt = 0; + if(!fTracks) + return pt; + if(matchIndex<0 || matchIndex>fTracks->GetEntries()){ + if(fDebug) + printf("track-matched index out of track array range!!!\n"); + return pt; + } + AliVTrack* track = static_cast(fTracks->At(matchIndex)); + if(!track){ + if(fDebug) + printf("track-matched pointer does not exist!!!\n"); + return pt; + } + if(fESD){ + if(!fPrTrCuts) + return pt; + if(!fPrTrCuts->IsSelected(track)) + return pt; + pt = track->Pt(); + } + return pt; +} +//________________________________________________________________________ +void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells() +{ + AliVCaloCells *cells = 0; + cells = fESDCells; + if (!cells) + cells = fAODCells; + if(!cells) + return; + Double_t maxe = 0; + Double_t maxphi = -10; + Int_t ncells = cells->GetNumberOfCells(); + Double_t eta,phi; + for (Int_t i=0; iGetCellNumber(i)); + Double_t e = cells->GetCellAmplitude(absid); + if(e>0.05) + fNCells50++; + else + continue; + fGeom->EtaPhiFromIndex(absid,eta,phi); + if(maxeFill(maxe,maxphi); + +} +//________________________________________________________________________ void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *) { // Called once at the end of the query.