| | 1 | // $Id$ |
| | 2 | |
| | 3 | #include "AliAnalysisTaskEMCALIsoPhoton.h" |
| | 4 | |
| | 5 | #include <TCanvas.h> |
| | 6 | #include <TChain.h> |
| | 7 | #include <TFile.h> |
| | 8 | #include <TH1F.h> |
| | 9 | #include <TH2F.h> |
| | 10 | #include <TH3F.h> |
| | 11 | #include <THnSparse.h> |
| | 12 | #include <TLorentzVector.h> |
| | 13 | #include <TList.h> |
| | 14 | |
| | 15 | #include "AliAnalysisManager.h" |
| | 16 | #include "AliAnalysisTask.h" |
| | 17 | #include "AliEMCALGeometry.h" |
| | 18 | #include "AliEMCALRecoUtils.h" |
| | 19 | #include "AliESDCaloCells.h" |
| | 20 | #include "AliESDCaloCluster.h" |
| | 21 | #include "AliESDEvent.h" |
| | 22 | #include "AliESDHeader.h" |
| | 23 | #include "AliESDInputHandler.h" |
| | 24 | #include "AliESDUtils.h" |
| | 25 | #include "AliESDtrack.h" |
| | 26 | #include "AliESDtrackCuts.h" |
| | 27 | #include "AliAODEvent.h" |
| | 28 | #include "AliAODTrack.h" |
| | 29 | #include "AliMCEvent.h" |
| | 30 | #include "AliMCEventHandler.h" |
| | 31 | #include "AliStack.h" |
| | 32 | #include "AliAODMCParticle.h" |
| | 33 | #include "AliVEvent.h" |
| | 34 | #include "AliVTrack.h" |
| | 35 | #include "AliV0vertexer.h" |
| | 36 | #include "AliVCluster.h" |
| | 37 | #include "AliOADBContainer.h" |
| | 38 | |
| | 39 | #include <iostream> |
| | 40 | using std::cout; |
| | 41 | using std::endl; |
| | 42 | |
| | 43 | ClassImp(AliAnalysisTaskEMCALIsoPhoton) |
| | 44 | |
| | 45 | //________________________________________________________________________ |
| | 46 | AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton() : |
| | 47 | AliAnalysisTaskSE(), |
| | 48 | fESDClusters(0), |
| | 49 | fAODClusters(0), |
| | 50 | fSelPrimTracks(0), |
| | 51 | fTracks(0), |
| | 52 | fAODMCParticles(0), |
| | 53 | fESDCells(0), |
| | 54 | fAODCells(0), |
| | 55 | fPrTrCuts(0), |
| | 56 | fCompTrCuts(0), |
| | 57 | fGeom(0x0), |
| | 58 | fGeoName("EMCAL_COMPLETEV1"), |
| | 59 | fOADBContainer(0), |
| | 60 | fVecPv(0.,0.,0.), |
| | 61 | fPeriod("LHC11c"), |
| | 62 | fTrigBit("kEMC7"), |
| | 63 | fIsTrain(0), |
| | 64 | fIsMc(0), |
| | 65 | fDebug(0), |
| | 66 | fPathStrOpt("/"), |
| | 67 | fExoticCut(0.97), |
| | 68 | fIsoConeR(0.4), |
| | 69 | fNDimensions(7), |
| | 70 | fECut(3.), |
| | 71 | fTrackMult(0), |
| | 72 | fMcIdFamily(""), |
| | 73 | fNClusForDirPho(0), |
| | 74 | fDirPhoPt(0), |
| | 75 | fHigherPtCone(0), |
| | 76 | fImportGeometryFromFile(0), |
| | 77 | fImportGeometryFilePath(""), |
| | 78 | fMaxPtTrack(0), |
| | 79 | fMaxEClus(0), |
| | 80 | fNCells50(0), |
| | 81 | fFilterBit(0), |
| | 82 | fSelHybrid(kFALSE), |
| | 83 | fFillQA(kFALSE), |
| | 84 | fClusIdFromTracks(""), |
| | 85 | fCpvFromTrack(kFALSE), |
| | 86 | fNBinsPt(200), |
| | 87 | fPtBinLowEdge(0), |
| | 88 | fPtBinHighEdge(100), |
| | 89 | fRemMatchClus(kFALSE), |
| | 90 | fMinIsoClusE(0), |
| | 91 | fNCuts(5), |
| | 92 | fTrCoreRem(kFALSE), |
| | 93 | fClusTDiff(30e-9), |
| | 94 | fPileUpRejSPD(kFALSE), |
| | 95 | fDistToBadChan(0), |
| | 96 | fInConeInvMass(""), |
| | 97 | fInConePairClEt(""), |
| | 98 | fESD(0), |
| | 99 | fAOD(0), |
| | 100 | fVEvent(0), |
| | 101 | fMCEvent(0), |
| | 102 | fStack(0), |
| | 103 | fOutputList(0), |
| | 104 | fEvtSel(0), |
| | 105 | fNClusEt10(0), |
| | 106 | fClusArrayNames(0), |
| | 107 | fRecoPV(0), |
| | 108 | fPVtxZ(0), |
| | 109 | fTrMultDist(0), |
| | 110 | fClusEtCPVSBGISO(0), |
| | 111 | fClusEtCPVBGISO(0), |
| | 112 | fMCDirPhotonPtEtaPhi(0), |
| | 113 | fMCIsoDirPhotonPtEtaPhi(0), |
| | 114 | fMCDirPhotonPtEtIso(0), |
| | 115 | fDecayPhotonPtMC(0), |
| | 116 | fCellAbsIdVsAmpl(0), |
| | 117 | fNClusHighClusE(0), |
| | 118 | fHigherPtConeM02(0), |
| | 119 | fClusEtMcPt(0), |
| | 120 | fClusMcDetaDphi(0), |
| | 121 | fNClusPerPho(0), |
| | 122 | fMcPtInConeBG(0), |
| | 123 | fMcPtInConeSBG(0), |
| | 124 | fMcPtInConeBGnoUE(0), |
| | 125 | fMcPtInConeSBGnoUE(0), |
| | 126 | fMcPtInConeTrBGnoUE(0), |
| | 127 | fMcPtInConeTrSBGnoUE(0), |
| | 128 | fMcPtInConeMcPhoPt(0), |
| | 129 | fAllIsoEtMcGamma(0), |
| | 130 | fAllIsoNoUeEtMcGamma(0), |
| | 131 | fMCDirPhotonPtEtaPhiNoClus(0), |
| | 132 | fEtCandIsoAndIsoWoPairEt(0), |
| | 133 | fInConePairedClusEtVsCandEt(0), |
| | 134 | fHnOutput(0), |
| | 135 | fQAList(0), |
| | 136 | fNTracks(0), |
| | 137 | fEmcNCells(0), |
| | 138 | fEmcNClus(0), |
| | 139 | fEmcNClusCut(0), |
| | 140 | fNTracksECut(0), |
| | 141 | fEmcNCellsCut(0), |
| | 142 | fEmcClusETM1(0), |
| | 143 | fEmcClusETM2(0), |
| | 144 | fEmcClusNotExo(0), |
| | 145 | fEmcClusEClusCuts(0), |
| | 146 | fEmcClusEPhi(0), |
| | 147 | fEmcClusEPhiCut(0), |
| | 148 | fEmcClusEEta(0), |
| | 149 | fEmcClusEEtaCut(0), |
| | 150 | fTrackPtPhi(0), |
| | 151 | fTrackPtPhiCut(0), |
| | 152 | fTrackPtEta(0), |
| | 153 | fTrackPtEtaCut(0), |
| | 154 | fMaxCellEPhi(0), |
| | 155 | fDetaDphiFromTM(0), |
| | 156 | fEoverPvsE(0) |
| | 157 | { |
| | 158 | // Default constructor. |
| | 159 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; |
| | 160 | } |
| | 161 | |
| | 162 | //________________________________________________________________________ |
| | 163 | AliAnalysisTaskEMCALIsoPhoton::AliAnalysisTaskEMCALIsoPhoton(const char *name) : |
| | 164 | AliAnalysisTaskSE(name), |
| | 165 | fESDClusters(0), |
| | 166 | fAODClusters(0), |
| | 167 | fSelPrimTracks(0), |
| | 168 | fTracks(0), |
| | 169 | fAODMCParticles(0), |
| | 170 | fESDCells(0), |
| | 171 | fAODCells(0), |
| | 172 | fPrTrCuts(0), |
| | 173 | fCompTrCuts(0), |
| | 174 | fGeom(0x0), |
| | 175 | fGeoName("EMCAL_COMPLETEV1"), |
| | 176 | fOADBContainer(0), |
| | 177 | fVecPv(0.,0.,0.), |
| | 178 | fPeriod("LHC11c"), |
| | 179 | fTrigBit("kEMC7"), |
| | 180 | fIsTrain(0), |
| | 181 | fIsMc(0), |
| | 182 | fDebug(0), |
| | 183 | fPathStrOpt("/"), |
| | 184 | fExoticCut(0.97), |
| | 185 | fIsoConeR(0.4), |
| | 186 | fNDimensions(7), |
| | 187 | fECut(3.), |
| | 188 | fTrackMult(0), |
| | 189 | fMcIdFamily(""), |
| | 190 | fNClusForDirPho(0), |
| | 191 | fDirPhoPt(0), |
| | 192 | fHigherPtCone(0), |
| | 193 | fImportGeometryFromFile(0), |
| | 194 | fImportGeometryFilePath(""), |
| | 195 | fMaxPtTrack(0), |
| | 196 | fMaxEClus(0), |
| | 197 | fNCells50(0), |
| | 198 | fFilterBit(0), |
| | 199 | fSelHybrid(kFALSE), |
| | 200 | fFillQA(kFALSE), |
| | 201 | fClusIdFromTracks(""), |
| | 202 | fCpvFromTrack(kFALSE), |
| | 203 | fNBinsPt(200), |
| | 204 | fPtBinLowEdge(0.), |
| | 205 | fPtBinHighEdge(100), |
| | 206 | fRemMatchClus(kFALSE), |
| | 207 | fMinIsoClusE(0), |
| | 208 | fNCuts(5), |
| | 209 | fTrCoreRem(kFALSE), |
| | 210 | fClusTDiff(30e-9), |
| | 211 | fPileUpRejSPD(kFALSE), |
| | 212 | fDistToBadChan(0), |
| | 213 | fInConeInvMass(""), |
| | 214 | fInConePairClEt(""), |
| | 215 | fESD(0), |
| | 216 | fAOD(0), |
| | 217 | fVEvent(0), |
| | 218 | fMCEvent(0), |
| | 219 | fStack(0), |
| | 220 | fOutputList(0), |
| | 221 | fEvtSel(0), |
| | 222 | fNClusEt10(0), |
| | 223 | fClusArrayNames(0), |
| | 224 | fRecoPV(0), |
| | 225 | fPVtxZ(0), |
| | 226 | fTrMultDist(0), |
| | 227 | fClusEtCPVSBGISO(0), |
| | 228 | fClusEtCPVBGISO(0), |
| | 229 | fMCDirPhotonPtEtaPhi(0), |
| | 230 | fMCIsoDirPhotonPtEtaPhi(0), |
| | 231 | fMCDirPhotonPtEtIso(0), |
| | 232 | fDecayPhotonPtMC(0), |
| | 233 | fCellAbsIdVsAmpl(0), |
| | 234 | fNClusHighClusE(0), |
| | 235 | fHigherPtConeM02(0), |
| | 236 | fClusEtMcPt(0), |
| | 237 | fClusMcDetaDphi(0), |
| | 238 | fNClusPerPho(0), |
| | 239 | fMcPtInConeBG(0), |
| | 240 | fMcPtInConeSBG(0), |
| | 241 | fMcPtInConeBGnoUE(0), |
| | 242 | fMcPtInConeSBGnoUE(0), |
| | 243 | fMcPtInConeTrBGnoUE(0), |
| | 244 | fMcPtInConeTrSBGnoUE(0), |
| | 245 | fMcPtInConeMcPhoPt(0), |
| | 246 | fAllIsoEtMcGamma(0), |
| | 247 | fAllIsoNoUeEtMcGamma(0), |
| | 248 | fMCDirPhotonPtEtaPhiNoClus(0), |
| | 249 | fEtCandIsoAndIsoWoPairEt(0), |
| | 250 | fInConePairedClusEtVsCandEt(0), |
| | 251 | fHnOutput(0), |
| | 252 | fQAList(0), |
| | 253 | fNTracks(0), |
| | 254 | fEmcNCells(0), |
| | 255 | fEmcNClus(0), |
| | 256 | fEmcNClusCut(0), |
| | 257 | fNTracksECut(0), |
| | 258 | fEmcNCellsCut(0), |
| | 259 | fEmcClusETM1(0), |
| | 260 | fEmcClusETM2(0), |
| | 261 | fEmcClusNotExo(0), |
| | 262 | fEmcClusEClusCuts(0), |
| | 263 | fEmcClusEPhi(0), |
| | 264 | fEmcClusEPhiCut(0), |
| | 265 | fEmcClusEEta(0), |
| | 266 | fEmcClusEEtaCut(0), |
| | 267 | fTrackPtPhi(0), |
| | 268 | fTrackPtPhiCut(0), |
| | 269 | fTrackPtEta(0), |
| | 270 | fTrackPtEtaCut(0), |
| | 271 | fMaxCellEPhi(0), |
| | 272 | fDetaDphiFromTM(0), |
| | 273 | fEoverPvsE(0) |
| | 274 | { |
| | 275 | // Constructor |
| | 276 | |
| | 277 | // Define input and output slots here |
| | 278 | // Input slot #0 works with a TChain |
| | 279 | DefineInput(0, TChain::Class()); |
| | 280 | // Output slot #0 id reserved by the base class for AOD |
| | 281 | // Output slot #1 writes into a TH1 container |
| | 282 | DefineOutput(1, TList::Class()); |
| | 283 | DefineOutput(2, TList::Class()); |
| | 284 | } |
| | 285 | |
| | 286 | //________________________________________________________________________ |
| | 287 | void AliAnalysisTaskEMCALIsoPhoton::UserCreateOutputObjects() |
| | 288 | { |
| | 289 | // Create histograms, called once. |
| | 290 | |
| | 291 | fESDClusters = new TObjArray(); |
| | 292 | fAODClusters = new TObjArray(); |
| | 293 | fSelPrimTracks = new TObjArray(); |
| | 294 | |
| | 295 | |
| | 296 | fOutputList = new TList(); |
| | 297 | fOutputList->SetOwner();// Container cleans up all histos (avoids leaks in merging) |
| | 298 | |
| | 299 | fGeom = AliEMCALGeometry::GetInstance(fGeoName.Data()); |
| | 300 | fOADBContainer = new AliOADBContainer("AliEMCALgeo"); |
| | 301 | fOADBContainer->InitFromFile(Form("$ALICE_ROOT/OADB/EMCAL/EMCALlocal2master.root"),"AliEMCALgeo"); |
| | 302 | |
| | 303 | fEvtSel = new TH1F("hEvtSel","Event selection counter (0=all trg, 1=pvz cut) ;evt cut ;dN/dcut}",2,0,2); |
| | 304 | fOutputList->Add(fEvtSel); |
| | 305 | |
| | 306 | fNClusEt10 = new TH1F("hNClusEt10","# of cluster with E_{T}>10 per event;E;",101,-0.5,100.5); |
| | 307 | fOutputList->Add(fNClusEt10); |
| | 308 | |
| | 309 | fClusArrayNames = new TH1F("hClusArrayNames","cluster array names (0=CaloClusters,1=EmcCaloClusters,2=Others);option;#events",3,0,3); |
| | 310 | fOutputList->Add(fClusArrayNames); |
| | 311 | |
| | 312 | fRecoPV = new TH1F("hRecoPV","Prim. vert. reconstruction (yes or no);reco (0=no, 1=yes) ;",2,-0.5,1.5); |
| | 313 | fOutputList->Add(fRecoPV); |
| | 314 | |
| | 315 | fPVtxZ = new TH1F("hPVtxZ","primary vertex Z before cut;prim-vz(cm) ;",200,-100,100); |
| | 316 | fOutputList->Add(fPVtxZ); |
| | 317 | |
| | 318 | fTrMultDist = new TH1F("fTrMultDist","track multiplicity;tracks/event;#events",200,0.5,200.5); |
| | 319 | fOutputList->Add(fTrMultDist); |
| | 320 | |
| | 321 | fClusEtCPVSBGISO = new TH2F("hClusEtCPVSBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.1<#lambda_{0}^{2}<0.3;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100); |
| | 322 | fClusEtCPVSBGISO->Sumw2(); |
| | 323 | fOutputList->Add(fClusEtCPVSBGISO); |
| | 324 | |
| | 325 | fClusEtCPVBGISO = new TH2F("hClusEtCPVBGISO","ISO^{TRK+EMC} vs. E_{T}^{clus} (after CPV and 0.5<#lambda_{0}^{2}<2.0;E_{T}^{clus} [GeV];ISO^{TRK+EMC} [GeV]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,100); |
| | 326 | fClusEtCPVBGISO->Sumw2(); |
| | 327 | fOutputList->Add(fClusEtCPVBGISO); |
| | 328 | |
| | 329 | fMCDirPhotonPtEtaPhi = new TH3F("hMCDirPhotonPtEtaPhi","photon (gq->#gammaq) p_{T}, #eta, #phi;GeV/c;#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20); |
| | 330 | fMCDirPhotonPtEtaPhi->Sumw2(); |
| | 331 | fOutputList->Add(fMCDirPhotonPtEtaPhi); |
| | 332 | |
| | 333 | fMCIsoDirPhotonPtEtaPhi = new TH3F("hMCIsoDirPhotonPtEtaPhi","photon (gq->#gammaq, isolated@MC) p_{T}, #eta, #phi;GeV/c;#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20); |
| | 334 | fMCIsoDirPhotonPtEtaPhi->Sumw2(); |
| | 335 | fOutputList->Add(fMCIsoDirPhotonPtEtaPhi); |
| | 336 | |
| | 337 | fMCDirPhotonPtEtIso = new TH2F("hMCDirPhotonPtEtIso",Form("photon (gq->#gammaq @MC) p_{T}, E_{T}^{ISO} (R=%1.1f);GeV/c;E_{T}^{ISO} GeV/c",fIsoConeR),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,20,-0.25,9.75); |
| | 338 | fMCDirPhotonPtEtIso->Sumw2(); |
| | 339 | fOutputList->Add(fMCDirPhotonPtEtIso); |
| | 340 | |
| | 341 | |
| | 342 | fDecayPhotonPtMC = new TH1F("hDecayPhotonPtMC","decay photon p_{T};GeV/c;dN/dp_{T} (c GeV^{-1})",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge); |
| | 343 | fDecayPhotonPtMC->Sumw2(); |
| | 344 | fOutputList->Add(fDecayPhotonPtMC); |
| | 345 | |
| | 346 | fCellAbsIdVsAmpl = new TH2F("hCellAbsIdVsAmpl","cell abs id vs cell amplitude (energy);E (GeV);ID",200,0,100,24*48*10,-0.5,24*48*10-0.5); |
| | 347 | fOutputList->Add(fCellAbsIdVsAmpl); |
| | 348 | |
| | 349 | fNClusHighClusE = new TH2F("hNClusHighClusE","total number of clusters vs. highest clus energy in the event;E (GeV);NClus",200,0,100,301,-0.5,300.5); |
| | 350 | fOutputList->Add(fNClusHighClusE); |
| | 351 | |
| | 352 | fHigherPtConeM02 = new TH2F("hHigherPtConeM02","#lambda_{0}^{2} vs. in-cone-p_{T}^{max};p_{T}^{max} (GeV/c, in the cone);#lambda_{0}^{2}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,400,0,4); |
| | 353 | fOutputList->Add(fHigherPtConeM02); |
| | 354 | |
| | 355 | fClusEtMcPt = new TH2F("hClusEtMcPt","E_{T}^{clus} vs. p_{T}^{mc}; p_{T}^{mc};E_{T}^{clus}",500,0,100,500,0,100); |
| | 356 | fOutputList->Add(fClusEtMcPt); |
| | 357 | |
| | 358 | fClusMcDetaDphi = new TH2F("hClusMcDetaDphi","#Delta#phi vs. #Delta#eta (reco-mc);#Delta#eta;#Delta#phi",100,-.7,.7,100,-.7,.7); |
| | 359 | fOutputList->Add(fClusMcDetaDphi); |
| | 360 | |
| | 361 | fNClusPerPho = new TH2F("hNClusPerPho","Number of clusters per prompt photon;p_{T}^{MC};N_{clus}",500,0,100,11,-0.5,10.5); |
| | 362 | fOutputList->Add(fNClusPerPho); |
| | 363 | |
| | 364 | fMcPtInConeBG = new TH2F("hMcPtInConeBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100); |
| | 365 | fOutputList->Add(fMcPtInConeBG); |
| | 366 | |
| | 367 | fMcPtInConeSBG = new TH2F("hMcPtInConeSBG","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100); |
| | 368 | fOutputList->Add(fMcPtInConeSBG); |
| | 369 | |
| | 370 | fMcPtInConeBGnoUE = new TH2F("hMcPtInConeBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (BG template);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100); |
| | 371 | fOutputList->Add(fMcPtInConeBGnoUE); |
| | 372 | |
| | 373 | fMcPtInConeSBGnoUE = new TH2F("hMcPtInConeSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK+EMC} (SBG range);ISO^{TRK+EMC} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100); |
| | 374 | fOutputList->Add(fMcPtInConeSBGnoUE); |
| | 375 | |
| | 376 | fMcPtInConeTrBGnoUE = new TH2F("hMcPtInConeTrBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (BG template);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100); |
| | 377 | fOutputList->Add(fMcPtInConeTrBGnoUE); |
| | 378 | |
| | 379 | fMcPtInConeTrSBGnoUE = new TH2F("hMcPtInConeTrSBGnoUE","#sum_{in-cone}p_{T}^{mc-primaries} vs. ISO^{TRK} (SBG range);ISO^{TRK} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",600,-10,50,1000,0,100); |
| | 380 | fOutputList->Add(fMcPtInConeTrSBGnoUE); |
| | 381 | |
| | 382 | fMcPtInConeMcPhoPt = new TH2F("hMcPtInConeMcPhoPt","#sum_{in-cone}p_{T}^{mc-primaries} vs. prompt photon p_{T};p_{T}^{mc-#gamma} (GeV);#sum_{in-cone}p_{T}^{mc-primaries}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,200,-0.25,99.75); |
| | 383 | fOutputList->Add(fMcPtInConeMcPhoPt); |
| | 384 | |
| | 385 | fAllIsoEtMcGamma = new TH2F("hAllIsoEtMcGammaE","ISO^{TRK+EMC} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50); |
| | 386 | fOutputList->Add(fAllIsoEtMcGamma); |
| | 387 | |
| | 388 | fAllIsoNoUeEtMcGamma = new TH2F("hAllIsoNoUeEtMcGammaE","ISO^{TRK+EMC}_{noue} vs. E_{T}^{clus} for clusters comming from a MC prompt #gamma; E_{T}^{clus} (GeV);ISO^{TRK+EMC}_{noue} (GeV);",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,600,-10,50); |
| | 389 | fOutputList->Add(fAllIsoNoUeEtMcGamma); |
| | 390 | |
| | 391 | |
| | 392 | fMCDirPhotonPtEtaPhiNoClus = new TH3F("hMCDirPhotonPhiEtaNoClus","p_{T}, #eta and #phi of prompt photons with no reco clusters;p_{T};#eta;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,154,-0.77,0.77,130,1.38,3.20); |
| | 393 | fOutputList->Add(fMCDirPhotonPtEtaPhiNoClus); |
| | 394 | |
| | 395 | fEtCandIsoAndIsoWoPairEt = new TH3F("hEtCandIsoAndIsoWoPairEt","E_{T}^{cand} vs. E_{T}^{ISO} (EMC+Trk) (0.1<M02<0.3, 0.110<m_{#gamma#gamma}<0.165 only);E_{T}^{cand}; E_{T}^{ISO}; E_{T}^{ISO} (w/o #pi^{0} pair E_{T})",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,1000,0,200,1000,0,200); |
| | 396 | fOutputList->Add(fEtCandIsoAndIsoWoPairEt); |
| | 397 | |
| | 398 | fInConePairedClusEtVsCandEt = new TH2F("hInConePairedClusEtVsCandEt","E_{T}^{partner} vs. E_{T}^{cand} (R<0.4, 0.110<m_{#gamma#gamma}<0.165);E_{T}^{cand};E_{T}^{partner}",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,200,0,40); |
| | 399 | fOutputList->Add(fInConePairedClusEtVsCandEt); |
| | 400 | |
| | 401 | Int_t nEt=fNBinsPt*5, 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=fNBinsPt, nInConeMass=100; |
| | 402 | Int_t bins[] = {nEt, nM02, nCeIso, nTrIso, nAllIso, nCeIsoNoUE, nAllIsoNoUE, nTrClDphi, nTrClDeta,nClEta,nClPhi,nTime,nMult,nPhoMcPt,nInConeMass}; |
| | 403 | fNDimensions = sizeof(bins)/sizeof(Int_t); |
| | 404 | const Int_t ndims = fNDimensions; |
| | 405 | Double_t xmin[] = { fPtBinLowEdge, 0., -10., -10., -10., -10., -10., -0.1,-0.05, -0.7, 1.4,-0.15e-06,-0.5,fPtBinLowEdge,0.0}; |
| | 406 | Double_t xmax[] = { fPtBinHighEdge, 4., 190., 190., 190., 190., 190., 0.1, 0.05, 0.7, 3.192, 0.15e-06,99.5,fPtBinHighEdge, 1.0}; |
| | 407 | if(fPeriod.Contains("11h")){ |
| | 408 | xmax[12]=3999.5; |
| | 409 | } |
| | 410 | 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); |
| | 411 | fHnOutput->Sumw2(); |
| | 412 | fOutputList->Add(fHnOutput); |
| | 413 | |
| | 414 | //QA outputs |
| | 415 | fQAList = new TList(); |
| | 416 | fQAList->SetOwner();// Container cleans up all histos (avoids leaks in merging) |
| | 417 | |
| | 418 | fNTracks = new TH1F("hNTracks","# of selected tracks;n-tracks;counts",120,-0.5,119.5); |
| | 419 | fNTracks->Sumw2(); |
| | 420 | fQAList->Add(fNTracks); |
| | 421 | |
| | 422 | fEmcNCells = new TH1F("fEmcNCells",";n/event;count",120,-0.5,119.5); |
| | 423 | fEmcNCells->Sumw2(); |
| | 424 | fQAList->Add(fEmcNCells); |
| | 425 | fEmcNClus = new TH1F("fEmcNClus",";n/event;count",120,-0.5,119.5); |
| | 426 | fEmcNClus->Sumw2(); |
| | 427 | fQAList->Add(fEmcNClus); |
| | 428 | fEmcNClusCut = new TH1F("fEmcNClusCut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5); |
| | 429 | fEmcNClusCut->Sumw2(); |
| | 430 | fQAList->Add(fEmcNClusCut); |
| | 431 | fNTracksECut = new TH1F("fNTracksECut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5); |
| | 432 | fNTracksECut->Sumw2(); |
| | 433 | fQAList->Add(fNTracksECut); |
| | 434 | fEmcNCellsCut = new TH1F("fEmcNCellsCut",Form("(at least one E_{clus}>%1.1f);n/event;count",fECut),120,-0.5,119.5); |
| | 435 | fEmcNCellsCut->Sumw2(); |
| | 436 | fQAList->Add(fEmcNCellsCut); |
| | 437 | fEmcClusETM1 = new TH1F("fEmcClusETM1","(method clus->GetTrackDx,z);GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge); |
| | 438 | fEmcClusETM1->Sumw2(); |
| | 439 | fQAList->Add(fEmcClusETM1); |
| | 440 | fEmcClusETM2 = new TH1F("fEmcClusETM2","(method track->GetEMCALcluster());GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge); |
| | 441 | fEmcClusETM2->Sumw2(); |
| | 442 | fQAList->Add(fEmcClusETM2); |
| | 443 | fEmcClusNotExo = new TH1F("fEmcClusNotExo","exotics removed;GeV;counts",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge); |
| | 444 | fEmcClusNotExo->Sumw2(); |
| | 445 | fQAList->Add(fEmcClusNotExo); |
| | 446 | fEmcClusEPhi = new TH2F("fEmcClusEPhi",";GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3); |
| | 447 | fEmcClusEPhi->Sumw2(); |
| | 448 | fQAList->Add(fEmcClusEPhi); |
| | 449 | fEmcClusEPhiCut = new TH2F("fEmcClusEPhiCut",Form("(at least one E_{clus}>%1.1f);GeV;#phi",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3); |
| | 450 | fEmcClusEPhiCut->Sumw2(); |
| | 451 | fQAList->Add(fEmcClusEPhiCut); |
| | 452 | fEmcClusEEta = new TH2F("fEmcClusEEta",";GeV;#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,19,-0.9,0.9); |
| | 453 | fEmcClusEEta->Sumw2(); |
| | 454 | fQAList->Add(fEmcClusEEta); |
| | 455 | fEmcClusEEtaCut = new TH2F("fEmcClusEEtaCut",Form("(at least one E_{clus}>%1.1f);GeV;#eta",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9); |
| | 456 | fEmcClusEEtaCut->Sumw2(); |
| | 457 | fQAList->Add(fEmcClusEEtaCut); |
| | 458 | fTrackPtPhi = new TH2F("fTrackPtPhi",";p_{T} [GeV/c];#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3); |
| | 459 | fTrackPtPhi->Sumw2(); |
| | 460 | fQAList->Add(fTrackPtPhi); |
| | 461 | fTrackPtPhiCut = new TH2F("fTrackPtPhiCut",Form("(at least one E_{clus}>%1.1f);p_{T} [GeV/c];#phi",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3); |
| | 462 | fTrackPtPhiCut->Sumw2(); |
| | 463 | fQAList->Add(fTrackPtPhiCut); |
| | 464 | fTrackPtEta = new TH2F("fTrackPtEta",";p_{T} [GeV/c];#eta",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9); |
| | 465 | fTrackPtEta->Sumw2(); |
| | 466 | fQAList->Add(fTrackPtEta); |
| | 467 | fTrackPtEtaCut = new TH2F("fTrackPtEtaCut",Form("(at least one E_{clus}>%1.1f);p_{T} [GeV/c];#eta",fECut),fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,18,-0.9,0.9); |
| | 468 | fTrackPtEtaCut->Sumw2(); |
| | 469 | fQAList->Add(fTrackPtEtaCut); |
| | 470 | fEmcClusEClusCuts = new TH2F("fEmcClusEClusCuts",";GeV;cut",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,fNCuts,-0.5,fNCuts-0.5); |
| | 471 | fEmcClusEClusCuts->Sumw2(); |
| | 472 | fQAList->Add(fEmcClusEClusCuts); |
| | 473 | |
| | 474 | fMaxCellEPhi = new TH2F("fMaxCellEPhi","Most energetic cell in event; GeV;#phi",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,63,0,6.3); |
| | 475 | fMaxCellEPhi->Sumw2(); |
| | 476 | fQAList->Add(fMaxCellEPhi); |
| | 477 | |
| | 478 | fDetaDphiFromTM = new TH2F("fDetaDphiFromTM","d#phi vs. d#eta of clusters from track->GetEMCALcluster();d#eta;d#phi",100,-0.05,0.05,200,-0.1,0.1); |
| | 479 | fDetaDphiFromTM->Sumw2(); |
| | 480 | fQAList->Add(fDetaDphiFromTM); |
| | 481 | |
| | 482 | fEoverPvsE = new TH2F("fEoverPvsE","E^{clus}/p^{track} vs E^{clus} (80<TPCsignal<100);E^{clus} [GeV];E^{clus}/p^{track} [c^{-1}]",fNBinsPt, fPtBinLowEdge,fPtBinHighEdge,100,0,2); |
| | 483 | fEoverPvsE->Sumw2(); |
| | 484 | fQAList->Add(fEoverPvsE); |
| | 485 | |
| | 486 | PostData(1, fOutputList); |
| | 487 | PostData(2, fQAList); |
| | 488 | } |
| | 489 | |
| | 490 | //________________________________________________________________________ |
| | 491 | void AliAnalysisTaskEMCALIsoPhoton::UserExec(Option_t *) |
| | 492 | { |
| | 493 | // Main loop, called for each event. |
| | 494 | fESDClusters = 0; |
| | 495 | fESDCells = 0; |
| | 496 | fAODClusters = 0; |
| | 497 | fAODCells = 0; |
| | 498 | // event trigger selection |
| | 499 | Bool_t isSelected = 0; |
| | 500 | if(fPeriod.Contains("11a")){ |
| | 501 | if(fTrigBit.Contains("kEMC")) |
| | 502 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC1); |
| | 503 | else |
| | 504 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); |
| | 505 | } |
| | 506 | else{ |
| | 507 | if(fTrigBit.Contains("kEMC")) |
| | 508 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kEMC7); |
| | 509 | else |
| | 510 | if(fTrigBit.Contains("kMB")) |
| | 511 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); |
| | 512 | else |
| | 513 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7); |
| | 514 | } |
| | 515 | if(fPeriod.Contains("11h")){ |
| | 516 | if(fTrigBit.Contains("kAny")) |
| | 517 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAny); |
| | 518 | if(fTrigBit.Contains("kAnyINT")) |
| | 519 | isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kAnyINT); |
| | 520 | } |
| | 521 | if(fIsMc) |
| | 522 | isSelected = kTRUE; |
| | 523 | if(fDebug) |
| | 524 | printf("isSelected = %d, fIsMC=%d\n", isSelected, fIsMc); |
| | 525 | if(!isSelected ) |
| | 526 | return; |
| | 527 | if(fIsMc){ |
| | 528 | TTree *tree = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetTree(); |
| | 529 | TFile *file = (TFile*)tree->GetCurrentFile(); |
| | 530 | TString filename = file->GetName(); |
| | 531 | if(!filename.Contains(fPathStrOpt.Data())) |
| | 532 | return; |
| | 533 | } |
| | 534 | fVEvent = (AliVEvent*)InputEvent(); |
| | 535 | if (!fVEvent) { |
| | 536 | printf("ERROR: event not available\n"); |
| | 537 | return; |
| | 538 | } |
| | 539 | Int_t runnumber = InputEvent()->GetRunNumber() ; |
| | 540 | if(fDebug) |
| | 541 | printf("run number = %d\n",runnumber); |
| | 542 | |
| | 543 | fESD = dynamic_cast<AliESDEvent*>(fVEvent); |
| | 544 | if(!fESD){ |
| | 545 | fAOD = dynamic_cast<AliAODEvent*>(fVEvent); |
| | 546 | if(!fAOD){ |
| | 547 | printf("ERROR: Invalid type of event!!!\n"); |
| | 548 | return; |
| | 549 | } |
| | 550 | else if(fDebug) |
| | 551 | printf("AOD EVENT!\n"); |
| | 552 | } |
| | 553 | |
| | 554 | fEvtSel->Fill(0); |
| | 555 | if(fDebug) |
| | 556 | printf("event is ok,\n run number=%d\n",runnumber); |
| | 557 | |
| | 558 | |
| | 559 | AliVVertex *pv = (AliVVertex*)fVEvent->GetPrimaryVertex(); |
| | 560 | Bool_t pvStatus = kTRUE; |
| | 561 | if(fESD){ |
| | 562 | AliESDVertex *esdv = (AliESDVertex*)fESD->GetPrimaryVertex(); |
| | 563 | pvStatus = esdv->GetStatus(); |
| | 564 | } |
| | 565 | /*if(fAOD){ |
| | 566 | AliAODVertex *aodv = (AliAODVertex*)fAOD->GetPrimaryVertex(); |
| | 567 | pvStatus = aodv->GetStatus(); |
| | 568 | }*/ |
| | 569 | if(!pv) |
| | 570 | return; |
| | 571 | if(!pvStatus) |
| | 572 | fRecoPV->Fill(0); |
| | 573 | else |
| | 574 | fRecoPV->Fill(1); |
| | 575 | fPVtxZ->Fill(pv->GetZ()); |
| | 576 | fVecPv.SetXYZ(pv->GetX(),pv->GetY(),pv->GetZ()); |
| | 577 | if(TMath::Abs(pv->GetZ())>10) |
| | 578 | return; |
| | 579 | if(fDebug) |
| | 580 | printf("passed vertex cut\n"); |
| | 581 | |
| | 582 | fEvtSel->Fill(1); |
| | 583 | if(fVEvent->IsPileupFromSPD(3, 0.8, 3., 2., 5.) && fPileUpRejSPD){ |
| | 584 | if(fDebug) |
| | 585 | printf("Event is SPD pile-up!***\n"); |
| | 586 | return; |
| | 587 | } |
| | 588 | if(fESD) |
| | 589 | fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("Tracks")); |
| | 590 | if(fAOD) |
| | 591 | fTracks = dynamic_cast<TClonesArray*>(fAOD->GetTracks()); |
| | 592 | |
| | 593 | if(!fTracks){ |
| | 594 | AliError("Track array in event is NULL!"); |
| | 595 | if(fDebug) |
| | 596 | printf("returning due to not finding tracks!\n"); |
| | 597 | return; |
| | 598 | } |
| | 599 | // Track loop to fill a pT spectrum |
| | 600 | const Int_t Ntracks = fTracks->GetEntriesFast(); |
| | 601 | for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) { |
| | 602 | // for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) { |
| | 603 | //AliESDtrack* track = (AliESDtrack*)fESD->GetTrack(iTracks); |
| | 604 | AliVTrack *track = (AliVTrack*)fTracks->At(iTracks); |
| | 605 | if (!track) |
| | 606 | continue; |
| | 607 | AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track); |
| | 608 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); |
| | 609 | if(esdTrack){ |
| | 610 | if(esdTrack->GetEMCALcluster()>0) |
| | 611 | fClusIdFromTracks.Append(Form("%d ",esdTrack->GetEMCALcluster())); |
| | 612 | if (fPrTrCuts && fPrTrCuts->IsSelected(track)){ |
| | 613 | fSelPrimTracks->Add(track); |
| | 614 | } else if(fCompTrCuts && fCompTrCuts->IsSelected(track)) { |
| | 615 | fSelPrimTracks->Add(track); |
| | 616 | } |
| | 617 | } |
| | 618 | else if(aodTrack){ |
| | 619 | if (fSelHybrid && !aodTrack->IsHybridGlobalConstrainedGlobal()) |
| | 620 | continue ; |
| | 621 | if(!fSelHybrid && !aodTrack->TestFilterBit(fFilterBit)) |
| | 622 | continue; |
| | 623 | fSelPrimTracks->Add(track); |
| | 624 | /*if(fTrackMaxPt<track->Pt()) |
| | 625 | fTrackMaxPt = track->Pt();*/ |
| | 626 | } |
| | 627 | } |
| | 628 | |
| | 629 | TObjArray *matEMCAL=(TObjArray*)fOADBContainer->GetObject(runnumber,"EmcalMatrices"); |
| | 630 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ |
| | 631 | if(fGeoName=="EMCAL_FIRSTYEARV1" && mod>3) |
| | 632 | break; |
| | 633 | /*if(fESD) |
| | 634 | fGeom->SetMisalMatrix(fESD->GetEMCALMatrix(mod), mod); |
| | 635 | else*/ |
| | 636 | // if(fVEvent->GetEMCALMatrix(mod)) |
| | 637 | fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod); |
| | 638 | fGeom->SetMisalMatrix(fGeomMatrix[mod] , mod); |
| | 639 | } |
| | 640 | |
| | 641 | if(fESD){ |
| | 642 | AliESDtrackCuts *fTrackCuts = new AliESDtrackCuts(); |
| | 643 | fTrackMult = fTrackCuts->GetReferenceMultiplicity(fESD);//kTrackletsITSTPC ,0.5); |
| | 644 | if(fDebug) |
| | 645 | printf("ESD Track mult= %d\n",fTrackMult); |
| | 646 | } |
| | 647 | else if(fAOD){ |
| | 648 | fTrackMult = Ntracks; |
| | 649 | if(fDebug) |
| | 650 | printf("AOD Track mult= %d\n",fTrackMult); |
| | 651 | } |
| | 652 | fTrMultDist->Fill(fTrackMult); |
| | 653 | TList *l = 0; |
| | 654 | TString clusArrayName = ""; |
| | 655 | if(fESD){ |
| | 656 | l = fESD->GetList(); |
| | 657 | if(fDebug) |
| | 658 | l->Print(); |
| | 659 | for(int nk=0;nk<l->GetEntries();nk++){ |
| | 660 | TObject *obj = (TObject*)l->At(nk); |
| | 661 | TString oname = obj->GetName(); |
| | 662 | if(oname.Contains("CaloClus")) |
| | 663 | clusArrayName = oname; |
| | 664 | else |
| | 665 | continue; |
| | 666 | if(clusArrayName=="CaloClusters") |
| | 667 | fClusArrayNames->Fill(0); |
| | 668 | else{ |
| | 669 | if(clusArrayName=="EmcCaloClusters") |
| | 670 | fClusArrayNames->Fill(1); |
| | 671 | else |
| | 672 | fClusArrayNames->Fill(2); |
| | 673 | } |
| | 674 | } |
| | 675 | fESDClusters = dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName)); |
| | 676 | fESDCells = fESD->GetEMCALCells(); |
| | 677 | if(fDebug) |
| | 678 | printf("ESD cluster mult= %d\n",fESDClusters->GetEntriesFast()); |
| | 679 | } |
| | 680 | else if(fAOD){ |
| | 681 | l = fAOD->GetList(); |
| | 682 | if(fDebug) |
| | 683 | l->Print(); |
| | 684 | //fAODClusters = dynamic_cast<TClonesArray*>(fAOD->GetCaloClusters()); |
| | 685 | for(int nk=0;nk<l->GetEntries();nk++){ |
| | 686 | TObject *obj = (TObject*)l->At(nk); |
| | 687 | TString oname = obj->GetName(); |
| | 688 | if(oname.Contains("aloClus")) |
| | 689 | clusArrayName = oname; |
| | 690 | else |
| | 691 | continue; |
| | 692 | if(clusArrayName=="caloClusters") |
| | 693 | fClusArrayNames->Fill(0); |
| | 694 | else{ |
| | 695 | if(clusArrayName=="EmcCaloClusters") |
| | 696 | fClusArrayNames->Fill(1); |
| | 697 | else |
| | 698 | fClusArrayNames->Fill(2); |
| | 699 | } |
| | 700 | } |
| | 701 | fAODClusters = dynamic_cast<TClonesArray*>(l->FindObject(clusArrayName)); |
| | 702 | fAODCells = fAOD->GetEMCALCells(); |
| | 703 | if(fDebug) |
| | 704 | printf("AOD cluster mult= %d\n",fAODClusters->GetEntriesFast()); |
| | 705 | } |
| | 706 | if(fDebug){ |
| | 707 | printf("clus array is named %s +++++++++\n",clusArrayName.Data()); |
| | 708 | } |
| | 709 | |
| | 710 | |
| | 711 | fMCEvent = MCEvent(); |
| | 712 | if(fMCEvent){ |
| | 713 | if(fDebug) |
| | 714 | std::cout<<"MCevent exists\n"; |
| | 715 | fStack = (AliStack*)fMCEvent->Stack(); |
| | 716 | if(!fStack) |
| | 717 | fAODMCParticles = (TClonesArray*)fVEvent->FindListObject("mcparticles"); |
| | 718 | } |
| | 719 | else{ |
| | 720 | if(fDebug) |
| | 721 | std::cout<<"ERROR: NO MC EVENT!!!!!!\n"; |
| | 722 | } |
| | 723 | LoopOnCells(); |
| | 724 | FollowGamma(); |
| | 725 | if(fDebug) |
| | 726 | printf("passed calling of FollowGamma\n"); |
| | 727 | FillClusHists(); |
| | 728 | if(fDebug) |
| | 729 | printf("passed calling of FillClusHists\n"); |
| | 730 | FillMcHists(); |
| | 731 | if(fDebug) |
| | 732 | printf("passed calling of FillMcHists\n"); |
| | 733 | if(fFillQA) |
| | 734 | FillQA(); |
| | 735 | if(fDebug) |
| | 736 | printf("passed calling of FillQA\n"); |
| | 737 | if(fESD) |
| | 738 | fESDClusters->Clear(); |
| | 739 | fSelPrimTracks->Clear(); |
| | 740 | fNClusForDirPho = 0; |
| | 741 | fNCells50 = 0; |
| | 742 | fClusIdFromTracks = ""; |
| | 743 | fVecPv.Clear(); |
| | 744 | |
| | 745 | PostData(1, fOutputList); |
| | 746 | PostData(2, fQAList); |
| | 747 | } |
| | 748 | |
| | 749 | //________________________________________________________________________ |
| | 750 | void AliAnalysisTaskEMCALIsoPhoton::FillClusHists() |
| | 751 | { |
| | 752 | if(fDebug) |
| | 753 | printf("Inside FillClusHists()....\n"); |
| | 754 | // Fill cluster histograms. |
| | 755 | TObjArray *clusters = fESDClusters; |
| | 756 | |
| | 757 | if (!clusters){ |
| | 758 | clusters = fAODClusters; |
| | 759 | if(fDebug) |
| | 760 | printf("ESD clusters empty..."); |
| | 761 | } |
| | 762 | if (!clusters){ |
| | 763 | if(fDebug) |
| | 764 | printf(" and AOD clusters as well!!!\n"); |
| | 765 | return; |
| | 766 | } |
| | 767 | if(fDebug) |
| | 768 | printf("\n"); |
| | 769 | |
| | 770 | const Int_t nclus = clusters->GetEntries(); |
| | 771 | if(nclus==0) |
| | 772 | return; |
| | 773 | if(fDebug) |
| | 774 | printf("Inside FillClusHists and there are %d clusters\n",nclus); |
| | 775 | Double_t maxE = 0; |
| | 776 | Int_t nclus10 = 0; |
| | 777 | Double_t ptmc=-1; |
| | 778 | for(Int_t ic=0;ic<nclus;ic++){ |
| | 779 | maxE=0; |
| | 780 | AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic)); |
| | 781 | if(!c){ |
| | 782 | if(fDebug) |
| | 783 | printf("cluster pointer does not exist! xxxx\n"); |
| | 784 | continue; |
| | 785 | } |
| | 786 | if(!c->IsEMCAL()){ |
| | 787 | if(fDebug) |
| | 788 | printf("cluster is not EMCAL! xxxx\n"); |
| | 789 | continue; |
| | 790 | } |
| | 791 | if(c->E()<fECut){ |
| | 792 | if(fDebug) |
| | 793 | printf("cluster has E<%1.1f! xxxx\n", fECut); |
| | 794 | continue; |
| | 795 | } |
| | 796 | if(fCpvFromTrack && fClusIdFromTracks.Contains(Form("%d",ic))){ |
| | 797 | if(fDebug) |
| | 798 | printf("cluster does not pass CPV criterion! xxxx\n"); |
| | 799 | continue; |
| | 800 | } |
| | 801 | if(IsExotic(c)){ |
| | 802 | if(fDebug) |
| | 803 | printf("cluster is exotic! xxxx\n"); |
| | 804 | continue; |
| | 805 | } |
| | 806 | if(c->GetDistanceToBadChannel()<fDistToBadChan){ |
| | 807 | if(fDebug) |
| | 808 | printf("cluster distance to bad channel is %1.1f (<%1.1f) xxxx\n",c->GetDistanceToBadChannel(),fDistToBadChan); |
| | 809 | continue; |
| | 810 | } |
| | 811 | Short_t id; |
| | 812 | Double_t Emax = GetMaxCellEnergy( c, id); |
| | 813 | if(fDebug) |
| | 814 | printf("cluster max cell E=%1.1f",Emax); |
| | 815 | Float_t clsPos[3] = {0,0,0}; |
| | 816 | c->GetPosition(clsPos); |
| | 817 | TVector3 clsVec(clsPos); |
| | 818 | clsVec -= fVecPv; |
| | 819 | Double_t Et = c->E()*TMath::Sin(clsVec.Theta()); |
| | 820 | if(fDebug) |
| | 821 | printf("\tcluster eta=%1.1f,phi=%1.1f,E=%1.1f\n",clsVec.Eta(),clsVec.Phi(),c->E()); |
| | 822 | if(Et>10) |
| | 823 | nclus10++; |
| | 824 | Float_t ceiso=0, cephiband=0, cecore=0; |
| | 825 | Float_t triso=0, trphiband=0, trcore=0; |
| | 826 | Float_t alliso=0, allphiband=0;//, allcore; |
| | 827 | Float_t phibandArea = (1.4 - 2*fIsoConeR)*2*fIsoConeR; |
| | 828 | Float_t netConeArea = TMath::Pi()*(fIsoConeR*fIsoConeR - 0.04*0.04); |
| | 829 | Bool_t isCPV = kFALSE; |
| | 830 | if(TMath::Abs(c->GetTrackDx())>0.03 || TMath::Abs(c->GetTrackDz())>0.02) |
| | 831 | isCPV = kTRUE; |
| | 832 | GetCeIso(clsVec, id, ceiso, cephiband, cecore, Et); |
| | 833 | GetTrIso(clsVec, triso, trphiband, trcore); |
| | 834 | Int_t nInConePairs = 0; |
| | 835 | Double_t onePairMass = 0; |
| | 836 | //--- |
| | 837 | //if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV){ |
| | 838 | TObjArray *inConeInvMassArr = (TObjArray*)fInConeInvMass.Tokenize(";"); |
| | 839 | TObjArray *inConePairClEt = (TObjArray*)fInConePairClEt.Tokenize(";"); |
| | 840 | nInConePairs = inConeInvMassArr->GetEntriesFast(); |
| | 841 | Int_t nInConePi0 = inConePairClEt->GetEntriesFast(); |
| | 842 | if(nInConePairs != nInConePi0) |
| | 843 | printf("Inconsistent number of in cone pairs!!!\n"); |
| | 844 | for(int ipair=0;ipair<nInConePairs;ipair++){ |
| | 845 | TObjString *obs = (TObjString*)inConeInvMassArr->At(ipair); |
| | 846 | TObjString *obet = (TObjString*)inConePairClEt->At(ipair); |
| | 847 | TString smass = obs->GetString(); |
| | 848 | TString spairEt = obet->GetString(); |
| | 849 | Double_t pairmass = smass.Atof(); |
| | 850 | Double_t pairEt = spairEt.Atof();//this must be zero when inv mass outside pi0 range |
| | 851 | if(0==ipair && nInConePairs==1) |
| | 852 | onePairMass = pairmass; |
| | 853 | if(fDebug) |
| | 854 | printf("=================+++++++++++++++Inv mass inside the cone for photon range: %1.1f,%1.1f,%1.1f+-++++-+-+-+-++-+-+-\n",Et,pairmass,ceiso+triso); |
| | 855 | fEtCandIsoAndIsoWoPairEt->Fill(Et,ceiso+triso,ceiso+triso-pairEt); |
| | 856 | } |
| | 857 | //} |
| | 858 | //--- |
| | 859 | Double_t dr = TMath::Sqrt(c->GetTrackDx()*c->GetTrackDx() + c->GetTrackDz()*c->GetTrackDz()); |
| | 860 | if(Et>10 && Et<15 && dr>0.025){ |
| | 861 | fHigherPtConeM02->Fill(fHigherPtCone,c->GetM02()); |
| | 862 | if(fDebug) |
| | 863 | printf("\t\tHigher track pt inside the cone: %1.1f GeV/c; M02=%1.2f\n",fHigherPtCone,c->GetM02()); |
| | 864 | } |
| | 865 | alliso = ceiso + triso; |
| | 866 | allphiband = cephiband + trphiband; |
| | 867 | //allcore = cecore + trcore; |
| | 868 | Float_t ceisoue = cephiband/phibandArea*netConeArea; |
| | 869 | Float_t trisoue = trphiband/phibandArea*netConeArea; |
| | 870 | Float_t allisoue = allphiband/phibandArea*netConeArea; |
| | 871 | Float_t mcptsum = GetMcPtSumInCone(clsVec.Eta(), clsVec.Phi(),fIsoConeR); |
| | 872 | if(fDebug && Et>10) |
| | 873 | printf("\t alliso=%1.1f, Et=%1.1f=-=-=-=-=\n",alliso,Et); |
| | 874 | if(c->GetM02()>0.5 && c->GetM02()<2.0){ |
| | 875 | fMcPtInConeBG->Fill(alliso-allisoue, mcptsum); |
| | 876 | fMcPtInConeBGnoUE->Fill(alliso, mcptsum); |
| | 877 | fMcPtInConeTrBGnoUE->Fill(triso, mcptsum); |
| | 878 | } |
| | 879 | if(c->GetM02()>0.1 && c->GetM02()<0.3 && dr>0.03){ |
| | 880 | fMcPtInConeSBG->Fill(alliso-allisoue, mcptsum); |
| | 881 | fMcPtInConeSBGnoUE->Fill(alliso, mcptsum); |
| | 882 | fMcPtInConeTrSBGnoUE->Fill(triso, mcptsum); |
| | 883 | if(fMcIdFamily.Contains((Form("%d",c->GetLabel())))){ |
| | 884 | fAllIsoEtMcGamma->Fill(Et, alliso-cecore-allisoue); |
| | 885 | fAllIsoNoUeEtMcGamma->Fill(Et, alliso-cecore); |
| | 886 | } |
| | 887 | } |
| | 888 | if(c->GetM02()>0.1 && c->GetM02()<0.3 && isCPV) |
| | 889 | fClusEtCPVSBGISO->Fill(Et,alliso - trcore); |
| | 890 | if(c->GetM02()>0.5 && c->GetM02()<2.0 && isCPV) |
| | 891 | fClusEtCPVBGISO->Fill(Et,alliso - trcore); |
| | 892 | const Int_t ndims = fNDimensions; |
| | 893 | Double_t outputValues[ndims]; |
| | 894 | if(mcptsum<2) |
| | 895 | ptmc = GetClusSource(c); |
| | 896 | else |
| | 897 | ptmc = -0.1; |
| | 898 | outputValues[0] = Et; |
| | 899 | outputValues[1] = c->GetM02(); |
| | 900 | outputValues[2] = ceiso/*cecore*/-ceisoue; |
| | 901 | outputValues[3] = triso-trisoue; |
| | 902 | outputValues[4] = alliso/*cecore*/-allisoue - trcore; |
| | 903 | outputValues[5] = ceiso; |
| | 904 | outputValues[6] = alliso - trcore; |
| | 905 | if(fDebug) |
| | 906 | printf("track-cluster dphi=%1.3f, deta=%1.3f\n",c->GetTrackDx(),c->GetTrackDz()); |
| | 907 | if(TMath::Abs(c->GetTrackDx())<0.1) |
| | 908 | outputValues[7] = c->GetTrackDx(); |
| | 909 | else |
| | 910 | outputValues[7] = 0.099*c->GetTrackDx()/TMath::Abs(c->GetTrackDx()); |
| | 911 | if(TMath::Abs(c->GetTrackDz())<0.05) |
| | 912 | outputValues[8] = c->GetTrackDz(); |
| | 913 | else |
| | 914 | outputValues[8] = 0.049*c->GetTrackDz()/TMath::Abs(c->GetTrackDz()); |
| | 915 | outputValues[9] = clsVec.Eta(); |
| | 916 | outputValues[10] = clsVec.Phi(); |
| | 917 | if(fESDCells) |
| | 918 | outputValues[11] = fESDCells->GetCellTime(id); |
| | 919 | else if(fAODCells) |
| | 920 | outputValues[11] = fAODCells->GetCellTime(id); |
| | 921 | outputValues[12] = fTrackMult; |
| | 922 | outputValues[13] = ptmc; |
| | 923 | if(nInConePairs == 1) |
| | 924 | outputValues[14] = onePairMass; |
| | 925 | else |
| | 926 | outputValues[14] = -1; |
| | 927 | fHnOutput->Fill(outputValues); |
| | 928 | if(c->E()>maxE) |
| | 929 | maxE = c->E(); |
| | 930 | } |
| | 931 | fNClusHighClusE->Fill(maxE,nclus); |
| | 932 | fMaxEClus = maxE; |
| | 933 | fNClusEt10->Fill(nclus10); |
| | 934 | fNClusPerPho->Fill(fDirPhoPt,fNClusForDirPho); |
| | 935 | } |
| | 936 | |
| | 937 | //________________________________________________________________________ |
| | 938 | void AliAnalysisTaskEMCALIsoPhoton::GetCeIso(TVector3 vec, Int_t maxid, Float_t &iso, Float_t &phiband, Float_t &core, Double_t EtCl) |
| | 939 | { |
| | 940 | if(fDebug) |
| | 941 | printf("....indside GetCeIso funtcion\n"); |
| | 942 | // Get cell isolation. |
| | 943 | AliVCaloCells *cells = fESDCells; |
| | 944 | if (!cells){ |
| | 945 | cells = fAODCells; |
| | 946 | if(fDebug) |
| | 947 | printf("ESD cells empty...\n"); |
| | 948 | } |
| | 949 | if (!cells){ |
| | 950 | if(fDebug) |
| | 951 | printf(" and AOD cells are empty as well!!!\n"); |
| | 952 | return; |
| | 953 | } |
| | 954 | |
| | 955 | TObjArray *clusters = fESDClusters; |
| | 956 | if (!clusters) |
| | 957 | clusters = fAODClusters; |
| | 958 | if (!clusters) |
| | 959 | return; |
| | 960 | |
| | 961 | fInConeInvMass = ""; |
| | 962 | fInConePairClEt=""; |
| | 963 | const Int_t nclus = clusters->GetEntries(); |
| | 964 | //const Int_t ncells = cells->GetNumberOfCells(); |
| | 965 | Float_t totiso=0; |
| | 966 | Float_t totband=0; |
| | 967 | Float_t totcore=0; |
| | 968 | Float_t etacl = vec.Eta(); |
| | 969 | Float_t phicl = vec.Phi(); |
| | 970 | if(phicl<0) |
| | 971 | phicl+=TMath::TwoPi(); |
| | 972 | /*Int_t absid = -1; |
| | 973 | Float_t eta=-1, phi=-1; |
| | 974 | for(int icell=0;icell<ncells;icell++){ |
| | 975 | absid = TMath::Abs(cells->GetCellNumber(icell)); |
| | 976 | Float_t celltime = cells->GetCellTime(absid); |
| | 977 | //if(TMath::Abs(celltime)>2e-8 && (!fIsMc)) |
| | 978 | if(TMath::Abs(celltime-maxtcl)>2e-8 ) |
| | 979 | continue; |
| | 980 | if(!fGeom) |
| | 981 | return; |
| | 982 | fGeom->EtaPhiFromIndex(absid,eta,phi); |
| | 983 | Float_t dphi = TMath::Abs(phi-phicl); |
| | 984 | Float_t deta = TMath::Abs(eta-etacl); |
| | 985 | Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);*/ |
| | 986 | for(int ic=0;ic<nclus;ic++){ |
| | 987 | AliVCluster *c = static_cast<AliVCluster*>(clusters->At(ic)); |
| | 988 | if(!c) |
| | 989 | continue; |
| | 990 | if(!c->IsEMCAL()) |
| | 991 | continue; |
| | 992 | if(c->E()<fMinIsoClusE) |
| | 993 | continue; |
| | 994 | Short_t id=-1; |
| | 995 | GetMaxCellEnergy( c, id); |
| | 996 | Double_t maxct = cells->GetCellTime(id); |
| | 997 | if(TMath::Abs(maxct)>fClusTDiff/*2.5e-9*/ && (!fIsMc)) |
| | 998 | continue; |
| | 999 | Float_t clsPos[3] = {0,0,0}; |
| | 1000 | c->GetPosition(clsPos); |
| | 1001 | TVector3 cv(clsPos); |
| | 1002 | cv -= fVecPv; |
| | 1003 | Double_t Et = c->E()*TMath::Sin(cv.Theta()); |
| | 1004 | Float_t dphi = (cv.Phi()-phicl); |
| | 1005 | Float_t deta = (cv.Eta()-etacl); |
| | 1006 | Float_t R = TMath::Sqrt(deta*deta + dphi*dphi); |
| | 1007 | if(R<0.007) |
| | 1008 | continue; |
| | 1009 | if(maxid==id) |
| | 1010 | continue; |
| | 1011 | Double_t matchedpt = GetTrackMatchedPt(c->GetTrackMatchedIndex()); |
| | 1012 | if(fCpvFromTrack){ |
| | 1013 | if(matchedpt>0 && fRemMatchClus ) |
| | 1014 | continue; |
| | 1015 | } else { |
| | 1016 | if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){ |
| | 1017 | if(fRemMatchClus){ |
| | 1018 | if(fDebug) |
| | 1019 | printf("This isolation cluster is matched to a track!++++++++++++++++++++++++++++++++++++++++++++++++++\n"); |
| | 1020 | continue; |
| | 1021 | } |
| | 1022 | } |
| | 1023 | } |
| | 1024 | Double_t nEt = TMath::Max(Et-matchedpt, 0.0); |
| | 1025 | if(nEt<0) |
| | 1026 | printf("nEt=%1.1f\n",nEt); |
| | 1027 | if(R<fIsoConeR){ |
| | 1028 | if(c->GetM02()>0.1 && c->GetM02()<0.3 && !(matchedpt>0)){ |
| | 1029 | TLorentzVector lv, lvec; |
| | 1030 | lv.SetPtEtaPhiM(Et,cv.Eta(),cv.Phi(),0); |
| | 1031 | lvec.SetPtEtaPhiM(EtCl,vec.Eta(),vec.Phi(),0); |
| | 1032 | TLorentzVector lpair = lv + lvec; |
| | 1033 | fInConeInvMass += Form("%f;",lpair.M()); |
| | 1034 | if(lpair.M()>0.11 && lpair.M()<0.165){ |
| | 1035 | fInConePairedClusEtVsCandEt->Fill(EtCl,Et); |
| | 1036 | fInConePairClEt += Form("%f;",Et); |
| | 1037 | //continue; |
| | 1038 | } |
| | 1039 | else |
| | 1040 | fInConePairClEt += Form("%f;",0.0); |
| | 1041 | /*if(lpair.M()>0.52 && lpair.M()<0.58) |
| | 1042 | continue;*/ |
| | 1043 | } |
| | 1044 | totiso += nEt; |
| | 1045 | if(R<0.04) |
| | 1046 | totcore += nEt; |
| | 1047 | } |
| | 1048 | else{ |
| | 1049 | if(dphi>fIsoConeR) |
| | 1050 | continue; |
| | 1051 | if(deta<fIsoConeR) |
| | 1052 | continue; |
| | 1053 | totband += nEt; |
| | 1054 | } |
| | 1055 | } |
| | 1056 | iso = totiso; |
| | 1057 | phiband = totband; |
| | 1058 | core = totcore; |
| | 1059 | } |
| | 1060 | //________________________________________________________________________ |
| | 1061 | void AliAnalysisTaskEMCALIsoPhoton::GetTrIso(TVector3 vec, Float_t &iso, Float_t &phiband, Float_t &core) |
| | 1062 | { |
| | 1063 | // Get track isolation. |
| | 1064 | |
| | 1065 | if(!fSelPrimTracks) |
| | 1066 | return; |
| | 1067 | fHigherPtCone = 0; |
| | 1068 | const Int_t ntracks = fSelPrimTracks->GetEntries(); |
| | 1069 | Double_t totiso=0; |
| | 1070 | Double_t totband=0; |
| | 1071 | Double_t totcore=0; |
| | 1072 | Double_t etacl = vec.Eta(); |
| | 1073 | Double_t phicl = vec.Phi(); |
| | 1074 | if(phicl<0) |
| | 1075 | phicl+=TMath::TwoPi(); |
| | 1076 | for(int itrack=0;itrack<ntracks;itrack++){ |
| | 1077 | AliVTrack *track = static_cast<AliVTrack*> (fSelPrimTracks->At(itrack)); |
| | 1078 | if(!track) |
| | 1079 | continue; |
| | 1080 | Double_t dphi = TMath::Abs(track->Phi()-phicl); |
| | 1081 | Double_t deta = TMath::Abs(track->Eta()-etacl); |
| | 1082 | Double_t R = TMath::Sqrt(deta*deta + dphi*dphi); |
| | 1083 | Double_t pt = track->Pt(); |
| | 1084 | if(pt>fHigherPtCone) |
| | 1085 | fHigherPtCone = pt; |
| | 1086 | if(R<fIsoConeR){ |
| | 1087 | totiso += track->Pt(); |
| | 1088 | if(R<0.04 && this->fTrCoreRem) |
| | 1089 | totcore += pt; |
| | 1090 | } |
| | 1091 | else{ |
| | 1092 | if(dphi>fIsoConeR) |
| | 1093 | continue; |
| | 1094 | if(deta<fIsoConeR) |
| | 1095 | continue; |
| | 1096 | totband += track->Pt(); |
| | 1097 | } |
| | 1098 | } |
| | 1099 | iso = totiso; |
| | 1100 | phiband = totband; |
| | 1101 | core = totcore; |
| | 1102 | } |
| | 1103 | |
| | 1104 | //________________________________________________________________________ |
| | 1105 | Double_t AliAnalysisTaskEMCALIsoPhoton::GetCrossEnergy(const AliVCluster *cluster, Short_t &idmax) |
| | 1106 | { |
| | 1107 | // Calculate the energy of cross cells around the leading cell. |
| | 1108 | |
| | 1109 | AliVCaloCells *cells = 0; |
| | 1110 | cells = fESDCells; |
| | 1111 | if (!cells) |
| | 1112 | cells = fAODCells; |
| | 1113 | if (!cells) |
| | 1114 | return 0; |
| | 1115 | |
| | 1116 | if (!fGeom) |
| | 1117 | return 0; |
| | 1118 | |
| | 1119 | Int_t iSupMod = -1; |
| | 1120 | Int_t iTower = -1; |
| | 1121 | Int_t iIphi = -1; |
| | 1122 | Int_t iIeta = -1; |
| | 1123 | Int_t iphi = -1; |
| | 1124 | Int_t ieta = -1; |
| | 1125 | Int_t iphis = -1; |
| | 1126 | Int_t ietas = -1; |
| | 1127 | |
| | 1128 | Double_t crossEnergy = 0; |
| | 1129 | |
| | 1130 | fGeom->GetCellIndex(idmax,iSupMod,iTower,iIphi,iIeta); |
| | 1131 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphis,ietas); |
| | 1132 | |
| | 1133 | Int_t ncells = cluster->GetNCells(); |
| | 1134 | for (Int_t i=0; i<ncells; i++) { |
| | 1135 | Int_t cellAbsId = cluster->GetCellAbsId(i); |
| | 1136 | fGeom->GetCellIndex(cellAbsId,iSupMod,iTower,iIphi,iIeta); |
| | 1137 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); |
| | 1138 | Int_t aphidiff = TMath::Abs(iphi-iphis); |
| | 1139 | if (aphidiff>1) |
| | 1140 | continue; |
| | 1141 | Int_t aetadiff = TMath::Abs(ieta-ietas); |
| | 1142 | if (aetadiff>1) |
| | 1143 | continue; |
| | 1144 | if ( (aphidiff==1 && aetadiff==0) || |
| | 1145 | (aphidiff==0 && aetadiff==1) ) { |
| | 1146 | crossEnergy += cells->GetCellAmplitude(cellAbsId); |
| | 1147 | } |
| | 1148 | } |
| | 1149 | |
| | 1150 | return crossEnergy; |
| | 1151 | } |
| | 1152 | |
| | 1153 | //________________________________________________________________________ |
| | 1154 | Double_t AliAnalysisTaskEMCALIsoPhoton ::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const |
| | 1155 | { |
| | 1156 | // Get maximum energy of attached cell. |
| | 1157 | |
| | 1158 | id = -1; |
| | 1159 | |
| | 1160 | AliVCaloCells *cells = 0; |
| | 1161 | cells = fESDCells; |
| | 1162 | if (!cells) |
| | 1163 | cells = fAODCells; |
| | 1164 | if(!cells) |
| | 1165 | return 0; |
| | 1166 | |
| | 1167 | Double_t maxe = 0; |
| | 1168 | Int_t ncells = cluster->GetNCells(); |
| | 1169 | for (Int_t i=0; i<ncells; i++) { |
| | 1170 | Double_t e = cells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); |
| | 1171 | if (e>maxe) { |
| | 1172 | maxe = e; |
| | 1173 | id = cluster->GetCellAbsId(i); |
| | 1174 | } |
| | 1175 | } |
| | 1176 | return maxe; |
| | 1177 | } |
| | 1178 | |
| | 1179 | //________________________________________________________________________ |
| | 1180 | void AliAnalysisTaskEMCALIsoPhoton ::FillMcHists() |
| | 1181 | { |
| | 1182 | if(!fStack && !fAODMCParticles) |
| | 1183 | return; |
| | 1184 | //ESD |
| | 1185 | if(fStack){ |
| | 1186 | Int_t nTracks = fStack->GetNtrack(); |
| | 1187 | if(fDebug) |
| | 1188 | printf("Inside FillMcHists and there are %d mcparts\n",nTracks); |
| | 1189 | for(Int_t iTrack=0;iTrack<nTracks;iTrack++){ |
| | 1190 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack)); |
| | 1191 | if(!mcp) |
| | 1192 | continue; |
| | 1193 | Int_t pdg = mcp->GetPdgCode(); |
| | 1194 | if(pdg!=22) |
| | 1195 | continue; |
| | 1196 | if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2) |
| | 1197 | continue; |
| | 1198 | Int_t imom = mcp->GetMother(0); |
| | 1199 | if(imom<0 || imom>nTracks) |
| | 1200 | continue; |
| | 1201 | TParticle *mcmom = static_cast<TParticle*>(fStack->Particle(imom)); |
| | 1202 | if(!mcmom) |
| | 1203 | continue; |
| | 1204 | Int_t pdgMom = mcmom->GetPdgCode(); |
| | 1205 | Double_t mcphi = mcp->Phi(); |
| | 1206 | Double_t mceta = mcp->Eta(); |
| | 1207 | if((imom==6 || imom==7) && pdgMom==22) { |
| | 1208 | fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
| | 1209 | Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR); |
| | 1210 | fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum); |
| | 1211 | if(ptsum<2) |
| | 1212 | fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
| | 1213 | if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR)) |
| | 1214 | fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum); |
| | 1215 | if(fNClusForDirPho==0) |
| | 1216 | fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
| | 1217 | if(fDebug){ |
| | 1218 | printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatusCode()); |
| | 1219 | printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatusCode()); |
| | 1220 | } |
| | 1221 | } |
| | 1222 | else{ |
| | 1223 | if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000) |
| | 1224 | fDecayPhotonPtMC->Fill(mcp->Pt()); |
| | 1225 | } |
| | 1226 | } |
| | 1227 | } |
| | 1228 | //AOD |
| | 1229 | else if(fAODMCParticles){ |
| | 1230 | Int_t nTracks = fAODMCParticles->GetEntriesFast(); |
| | 1231 | if(fDebug) |
| | 1232 | printf("Inside FillMcHists and there are %d mcparts\n",nTracks); |
| | 1233 | for(Int_t iTrack=0;iTrack<nTracks;iTrack++){ |
| | 1234 | AliAODMCParticle *mcp = dynamic_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack)); |
| | 1235 | if(!mcp) |
| | 1236 | continue; |
| | 1237 | Int_t pdg = mcp->GetPdgCode(); |
| | 1238 | if(pdg!=22) |
| | 1239 | continue; |
| | 1240 | if(TMath::Abs(mcp->Eta())>0.7 ||mcp->Phi()<1.4 || mcp->Phi()>3.2) |
| | 1241 | continue; |
| | 1242 | Int_t imom = mcp->GetMother(); |
| | 1243 | if(imom<0 || imom>nTracks) |
| | 1244 | continue; |
| | 1245 | AliAODMCParticle *mcmom = static_cast<AliAODMCParticle*>(fAODMCParticles->At(imom)); |
| | 1246 | if(!mcmom) |
| | 1247 | continue; |
| | 1248 | Int_t pdgMom = mcmom->GetPdgCode(); |
| | 1249 | Double_t mcphi = mcp->Phi(); |
| | 1250 | Double_t mceta = mcp->Eta(); |
| | 1251 | if((imom==6 || imom==7) && pdgMom==22) { |
| | 1252 | fMCDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
| | 1253 | Float_t ptsum = GetMcPtSumInCone(mcp->Eta(), mcp->Phi(), fIsoConeR); |
| | 1254 | fMcPtInConeMcPhoPt->Fill(mcp->Pt(),ptsum); |
| | 1255 | if(ptsum<2) |
| | 1256 | fMCIsoDirPhotonPtEtaPhi->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
| | 1257 | if(mcphi<(3.14-fIsoConeR) && mcphi>(1.4+fIsoConeR) && TMath::Abs(mceta)<(0.7-fIsoConeR)) |
| | 1258 | fMCDirPhotonPtEtIso->Fill(mcp->Pt(),ptsum); |
| | 1259 | if(fNClusForDirPho==0) |
| | 1260 | fMCDirPhotonPtEtaPhiNoClus->Fill(mcp->Pt(),mcp->Eta(),mcp->Phi()); |
| | 1261 | if(fDebug){ |
| | 1262 | printf("Found \"photonic\" parton at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d,\n",imom,mcmom->Pt(), mcmom->Eta(), mcmom->Phi(), mcmom->GetStatus()); |
| | 1263 | printf("with a final photon at position %d, with pt=%1.1f, eta=%1.1f and phi=%1.1f, and status=%d\n",iTrack,mcp->Pt(), mcp->Eta(), mcp->Phi(),mcp->GetStatus()); |
| | 1264 | } |
| | 1265 | } |
| | 1266 | else{ |
| | 1267 | if(TMath::Abs(pdgMom)>100 && TMath::Abs(pdgMom)<1000) |
| | 1268 | fDecayPhotonPtMC->Fill(mcp->Pt()); |
| | 1269 | } |
| | 1270 | } |
| | 1271 | } |
| | 1272 | } |
| | 1273 | //________________________________________________________________________ |
| | 1274 | Float_t AliAnalysisTaskEMCALIsoPhoton::GetClusSource(const AliVCluster *c) |
| | 1275 | { |
| | 1276 | if(!c) |
| | 1277 | return -0.1; |
| | 1278 | if(!fStack && !fAODMCParticles) |
| | 1279 | return -0.1; |
| | 1280 | Int_t clabel = c->GetLabel(); |
| | 1281 | if(fDebug && fMcIdFamily.Contains(Form("%d",clabel))) |
| | 1282 | printf("\n\t ==== Label %d is a descendent of the prompt photon ====\n\n",clabel); |
| | 1283 | if(!fMcIdFamily.Contains(Form("%d",clabel))) |
| | 1284 | return -0.1; |
| | 1285 | fNClusForDirPho++; |
| | 1286 | TString partonposstr = (TSubString)fMcIdFamily.operator()(0,1); |
| | 1287 | Int_t partonpos = partonposstr.Atoi(); |
| | 1288 | if(fDebug) |
| | 1289 | printf("\t^^^^ parton position = %d ^^^^\n",partonpos); |
| | 1290 | Float_t clsPos[3] = {0,0,0}; |
| | 1291 | c->GetPosition(clsPos); |
| | 1292 | TVector3 clsVec(clsPos); |
| | 1293 | clsVec -= fVecPv; |
| | 1294 | Double_t Et = c->E()*TMath::Sin(clsVec.Theta()); |
| | 1295 | //ESD |
| | 1296 | if(fStack){ |
| | 1297 | Int_t nmcp = fStack->GetNtrack(); |
| | 1298 | if(clabel<0 || clabel>nmcp) |
| | 1299 | return -0.1; |
| | 1300 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(partonpos)); |
| | 1301 | if(!mcp) |
| | 1302 | return -0.1; |
| | 1303 | if(fDebug){ |
| | 1304 | 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); |
| | 1305 | } |
| | 1306 | Int_t lab1 = mcp->GetFirstDaughter(); |
| | 1307 | if(lab1<0 || lab1>nmcp) |
| | 1308 | return -0.1; |
| | 1309 | TParticle *mcd = static_cast<TParticle*>(fStack->Particle(lab1)); |
| | 1310 | if(!mcd) |
| | 1311 | return -0.1; |
| | 1312 | if(fDebug) |
| | 1313 | 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); |
| | 1314 | if(mcd->GetPdgCode()==22){ |
| | 1315 | fClusEtMcPt->Fill(mcd->Pt(), Et); |
| | 1316 | fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi()); |
| | 1317 | } |
| | 1318 | else{ |
| | 1319 | if(fDebug) |
| | 1320 | printf("Warning: daughter of photon parton is not a photon\n"); |
| | 1321 | return -0.1; |
| | 1322 | } |
| | 1323 | } |
| | 1324 | //AOD |
| | 1325 | else if(fAODMCParticles){ |
| | 1326 | Int_t nmcp = fAODMCParticles->GetEntriesFast(); |
| | 1327 | if(clabel<0 || clabel>nmcp) |
| | 1328 | return -0.1; |
| | 1329 | AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(partonpos)); |
| | 1330 | if(!mcp) |
| | 1331 | return -0.1; |
| | 1332 | if(fDebug){ |
| | 1333 | printf("\tclus mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),clabel); |
| | 1334 | } |
| | 1335 | Int_t lab1 = mcp->GetDaughter(0); |
| | 1336 | if(lab1<0 || lab1>nmcp) |
| | 1337 | return -0.1; |
| | 1338 | AliAODMCParticle *mcd = static_cast<AliAODMCParticle*>(fAODMCParticles->At(lab1)); |
| | 1339 | if(!mcd) |
| | 1340 | return -0.1; |
| | 1341 | if(fDebug) |
| | 1342 | printf("\t\tmom mc truth eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, stackpos=%d\n",mcd->Eta(),mcd->Phi(),mcd->E(),mcd->GetPdgCode(),lab1); |
| | 1343 | if(mcd->GetPdgCode()==22){ |
| | 1344 | fClusEtMcPt->Fill(mcd->Pt(), Et); |
| | 1345 | fClusMcDetaDphi->Fill(clsVec.Eta() - mcd->Eta(), clsVec.Phi() - mcd->Phi()); |
| | 1346 | } |
| | 1347 | else{ |
| | 1348 | printf("Warning: daughter of photon parton is not a photon\n"); |
| | 1349 | return -0.1; |
| | 1350 | } |
| | 1351 | } |
| | 1352 | return fDirPhoPt; |
| | 1353 | } |
| | 1354 | //________________________________________________________________________ |
| | 1355 | void AliAnalysisTaskEMCALIsoPhoton::FollowGamma() |
| | 1356 | { |
| | 1357 | if(!fStack && !fAODMCParticles) |
| | 1358 | return; |
| | 1359 | Int_t selfid = 6; |
| | 1360 | //ESD |
| | 1361 | if(fStack){ |
| | 1362 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(selfid)); |
| | 1363 | if(!mcp) |
| | 1364 | return; |
| | 1365 | if(mcp->GetPdgCode()!=22){ |
| | 1366 | mcp = static_cast<TParticle*>(fStack->Particle(++selfid)); |
| | 1367 | if(!mcp) |
| | 1368 | return; |
| | 1369 | } |
| | 1370 | Int_t daug0f = mcp->GetFirstDaughter(); |
| | 1371 | Int_t daug0l = mcp->GetLastDaughter(); |
| | 1372 | Int_t nd0 = daug0l - daug0f; |
| | 1373 | if(fDebug) |
| | 1374 | printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->Energy(),mcp->GetPdgCode(),nd0+1); |
| | 1375 | fMcIdFamily = Form("%d,",selfid); |
| | 1376 | GetDaughtersInfo(daug0f,daug0l, selfid,""); |
| | 1377 | if(fDebug){ |
| | 1378 | printf("\t---- end of the prompt gamma's family tree ----\n\n"); |
| | 1379 | printf("Family id string = %s,\n\n",fMcIdFamily.Data()); |
| | 1380 | } |
| | 1381 | TParticle *md = static_cast<TParticle*>(fStack->Particle(daug0f)); |
| | 1382 | if(!md) |
| | 1383 | return; |
| | 1384 | fDirPhoPt = md->Pt(); |
| | 1385 | } |
| | 1386 | //AOD |
| | 1387 | else if(fAODMCParticles){ |
| | 1388 | AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(selfid)); |
| | 1389 | if(!mcp) |
| | 1390 | return; |
| | 1391 | if(mcp->GetPdgCode()!=22){ |
| | 1392 | mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(++selfid)); |
| | 1393 | if(!mcp) |
| | 1394 | return; |
| | 1395 | } |
| | 1396 | Int_t daug0f = mcp->GetDaughter(0); |
| | 1397 | Int_t daug0l = mcp->GetDaughter(mcp->GetNDaughters()-1); |
| | 1398 | Int_t nd0 = daug0l - daug0f; |
| | 1399 | if(fDebug) |
| | 1400 | printf("\n\tGenerated gamma (%d) eta=%1.1f,phi=%1.1f,E=%1.1f, pdgcode=%d, n-daug=%d\n",selfid,mcp->Eta(),mcp->Phi(),mcp->E(),mcp->GetPdgCode(),nd0+1); |
| | 1401 | fMcIdFamily = Form("%d,",selfid); |
| | 1402 | GetDaughtersInfo(daug0f,daug0l, selfid,""); |
| | 1403 | if(fDebug){ |
| | 1404 | printf("\t---- end of the prompt gamma's family tree ----\n\n"); |
| | 1405 | printf("Family id string = %s,\n\n",fMcIdFamily.Data()); |
| | 1406 | } |
| | 1407 | AliAODMCParticle *md = static_cast<AliAODMCParticle*>(fAODMCParticles->At(daug0f)); |
| | 1408 | if(!md) |
| | 1409 | return; |
| | 1410 | fDirPhoPt = md->Pt(); |
| | 1411 | } |
| | 1412 | |
| | 1413 | } |
| | 1414 | //________________________________________________________________________ |
| | 1415 | void AliAnalysisTaskEMCALIsoPhoton::GetDaughtersInfo(int firstd, int lastd, int selfid, const char *inputind) |
| | 1416 | { |
| | 1417 | if(fStack){ |
| | 1418 | int nmcp = fStack->GetNtrack(); |
| | 1419 | if(firstd<0 || firstd>nmcp) |
| | 1420 | return; |
| | 1421 | if(lastd<0 || lastd>nmcp) |
| | 1422 | return; |
| | 1423 | if(firstd>lastd){ |
| | 1424 | printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd); |
| | 1425 | return; |
| | 1426 | } |
| | 1427 | TString indenter = Form("\t%s",inputind); |
| | 1428 | TParticle *dp = 0x0; |
| | 1429 | if(fDebug) |
| | 1430 | printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid); |
| | 1431 | for(int id=firstd; id<lastd+1; id++){ |
| | 1432 | dp = static_cast<TParticle*>(fStack->Particle(id)); |
| | 1433 | if(!dp) |
| | 1434 | continue; |
| | 1435 | Int_t dfd = dp->GetFirstDaughter(); |
| | 1436 | Int_t dld = dp->GetLastDaughter(); |
| | 1437 | Int_t dnd = dld - dfd + 1; |
| | 1438 | Float_t vr = TMath::Sqrt(dp->Vx()*dp->Vx()+dp->Vy()*dp->Vy()); |
| | 1439 | if(fDebug) |
| | 1440 | printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->Energy(), vr, dp->GetPdgCode(), dnd, dfd, dld); |
| | 1441 | fMcIdFamily += Form("%d,",id); |
| | 1442 | GetDaughtersInfo(dfd,dld,id,indenter.Data()); |
| | 1443 | } |
| | 1444 | } |
| | 1445 | if(fAODMCParticles){ |
| | 1446 | int nmcp = fAODMCParticles->GetEntriesFast(); |
| | 1447 | if(firstd<0 || firstd>nmcp) |
| | 1448 | return; |
| | 1449 | if(lastd<0 || lastd>nmcp) |
| | 1450 | return; |
| | 1451 | if(firstd>lastd){ |
| | 1452 | printf("WARNING: First daughter > last (%d,%d)\n",firstd,lastd); |
| | 1453 | return; |
| | 1454 | } |
| | 1455 | TString indenter = Form("\t%s",inputind); |
| | 1456 | AliAODMCParticle *dp = 0x0; |
| | 1457 | if(fDebug) |
| | 1458 | printf("\t%s--- Daughters of particle %d ---\n", indenter.Data(), selfid); |
| | 1459 | for(int id=firstd; id<lastd+1; id++){ |
| | 1460 | dp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(id)); |
| | 1461 | if(!dp) |
| | 1462 | continue; |
| | 1463 | Int_t dfd = dp->GetDaughter(0); |
| | 1464 | Int_t dld = dp->GetDaughter(dp->GetNDaughters()-1); |
| | 1465 | Int_t dnd = dld - dfd + 1; |
| | 1466 | Float_t vr = TMath::Sqrt(dp->Xv()*dp->Xv()+dp->Xv()*dp->Xv()); |
| | 1467 | if(fDebug) |
| | 1468 | printf("\t%sParticle daughter(%d) eta=%1.1f,phi=%1.1f,E=%1.1f, vR=%1.1f, pdgcode=%d, n-daug=%d(%d,%d)\n", indenter.Data(),id, dp->Eta(), dp->Phi(), dp->E(), vr, dp->GetPdgCode(), dnd, dfd, dld); |
| | 1469 | fMcIdFamily += Form("%d,",id); |
| | 1470 | GetDaughtersInfo(dfd,dld,id,indenter.Data()); |
| | 1471 | } |
| | 1472 | } |
| | 1473 | } |
| | 1474 | |
| | 1475 | //________________________________________________________________________ |
| | 1476 | Float_t AliAnalysisTaskEMCALIsoPhoton::GetMcPtSumInCone(Float_t etaclus, Float_t phiclus, Float_t R) |
| | 1477 | { |
| | 1478 | if(!fStack && !fAODMCParticles) |
| | 1479 | return 0; |
| | 1480 | if(fDebug) |
| | 1481 | printf("Inside GetMcPtSumInCone!!\n"); |
| | 1482 | Float_t ptsum = 0; |
| | 1483 | TString addpartlabels = ""; |
| | 1484 | //ESD |
| | 1485 | if(fStack){ |
| | 1486 | Int_t nTracks = fStack->GetNtrack(); |
| | 1487 | for(Int_t iTrack=9;iTrack<nTracks;iTrack++){ |
| | 1488 | TParticle *mcp = static_cast<TParticle*>(fStack->Particle(iTrack)); |
| | 1489 | if(!mcp) |
| | 1490 | continue; |
| | 1491 | Int_t status = mcp->GetStatusCode(); |
| | 1492 | if(status!=1) |
| | 1493 | continue; |
| | 1494 | Float_t mcvr = TMath::Sqrt(mcp->Vx()*mcp->Vx()+ mcp->Vy()* mcp->Vy() + mcp->Vz()*mcp->Vz()); |
| | 1495 | if(mcvr>1) |
| | 1496 | continue; |
| | 1497 | /*else { |
| | 1498 | if(fDebug) |
| | 1499 | printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz()); |
| | 1500 | }*/ |
| | 1501 | Float_t dphi = mcp->Phi() - phiclus; |
| | 1502 | Float_t deta = mcp->Eta() - etaclus; |
| | 1503 | if(fDebug && TMath::Abs(dphi)<0.01) |
| | 1504 | printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta); |
| | 1505 | |
| | 1506 | if(deta>R || dphi>R) |
| | 1507 | continue; |
| | 1508 | Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta); |
| | 1509 | if(dR>R) |
| | 1510 | continue; |
| | 1511 | addpartlabels += Form("%d,",iTrack); |
| | 1512 | if(mcp->Pt()<0.2) |
| | 1513 | continue; |
| | 1514 | ptsum += mcp->Pt(); |
| | 1515 | } |
| | 1516 | } |
| | 1517 | //AOD |
| | 1518 | if(fAODMCParticles){ |
| | 1519 | Int_t nTracks = fAODMCParticles->GetEntriesFast(); |
| | 1520 | for(Int_t iTrack=9;iTrack<nTracks;iTrack++){ |
| | 1521 | AliAODMCParticle *mcp = static_cast<AliAODMCParticle*>(fAODMCParticles->At(iTrack)); |
| | 1522 | if(!mcp) |
| | 1523 | continue; |
| | 1524 | Int_t status = mcp->GetStatus(); |
| | 1525 | if(status!=1) |
| | 1526 | continue; |
| | 1527 | Float_t mcvr = TMath::Sqrt(mcp->Xv()*mcp->Xv()+ mcp->Yv()* mcp->Yv() + mcp->Zv()*mcp->Zv()); |
| | 1528 | if(mcvr>1) |
| | 1529 | continue; |
| | 1530 | /*else { |
| | 1531 | if(fDebug) |
| | 1532 | printf(" >>>> mcp Rho, Vx, Vy, Vz = %1.1f,%1.1f,%1.1f,%1.1f.......\n",mcp->Rho(),mcp->Vx(), mcp->Vy(),mcp->Vz()); |
| | 1533 | }*/ |
| | 1534 | Float_t dphi = mcp->Phi() - phiclus; |
| | 1535 | Float_t deta = mcp->Eta() - etaclus; |
| | 1536 | if(fDebug && TMath::Abs(dphi)<0.01) |
| | 1537 | printf(" >>>> mcphi = %1.1f, mceta = %1.1f\n>>>> dphi = %1.1f, deta = %1.1f\n", mcp->Phi(), mcp->Eta(),dphi,deta); |
| | 1538 | |
| | 1539 | if(deta>R || dphi>R) |
| | 1540 | continue; |
| | 1541 | Float_t dR = TMath::Sqrt(dphi*dphi + deta*deta); |
| | 1542 | if(dR>R) |
| | 1543 | continue; |
| | 1544 | addpartlabels += Form("%d,",iTrack); |
| | 1545 | if(mcp->Pt()<0.2) |
| | 1546 | continue; |
| | 1547 | ptsum += mcp->Pt(); |
| | 1548 | } |
| | 1549 | } |
| | 1550 | return ptsum; |
| | 1551 | } |
| | 1552 | //________________________________________________________________________ |
| | 1553 | void AliAnalysisTaskEMCALIsoPhoton::FillQA() |
| | 1554 | { |
| | 1555 | |
| | 1556 | TObjArray *clusters = fESDClusters; |
| | 1557 | //"none", "exotic", "exo+cpv1", "exo+cpv1+time", "exo+cpv1+time+m02"), |
| | 1558 | if (!clusters){ |
| | 1559 | clusters = fAODClusters; |
| | 1560 | if(fDebug) |
| | 1561 | printf("ESD clusters empty..."); |
| | 1562 | } |
| | 1563 | if (!clusters){ |
| | 1564 | if(fDebug) |
| | 1565 | printf(" and AOD clusters as well!!!\n"); |
| | 1566 | return; |
| | 1567 | } |
| | 1568 | if(!fSelPrimTracks) |
| | 1569 | return; |
| | 1570 | const int ntracks = fSelPrimTracks->GetEntriesFast(); |
| | 1571 | const int ncells = fNCells50;//fESDCells->GetNumberOfCells(); |
| | 1572 | const Int_t nclus = clusters->GetEntries(); |
| | 1573 | fNTracks->Fill(ntracks); |
| | 1574 | fEmcNCells->Fill(ncells); |
| | 1575 | fEmcNClus->Fill(nclus); |
| | 1576 | if(fMaxEClus>fECut){ |
| | 1577 | fNTracksECut->Fill(ntracks); |
| | 1578 | fEmcNCellsCut->Fill(ncells); |
| | 1579 | fEmcNClusCut->Fill(nclus); |
| | 1580 | } |
| | 1581 | for(int it=0;it<ntracks;it++){ |
| | 1582 | AliVTrack *t = (AliVTrack*)fSelPrimTracks->At(it); |
| | 1583 | if(!t) |
| | 1584 | continue; |
| | 1585 | fTrackPtPhi->Fill(t->Pt(),t->Phi()); |
| | 1586 | fTrackPtEta->Fill(t->Pt(),t->Eta()); |
| | 1587 | if(fMaxEClus>fECut){ |
| | 1588 | fTrackPtPhiCut->Fill(t->Pt(), t->Phi()); |
| | 1589 | fTrackPtEtaCut->Fill(t->Pt(), t->Eta()); |
| | 1590 | } |
| | 1591 | if(t->GetTPCsignal()<80 || t->GetTPCsignal()>100) |
| | 1592 | continue; |
| | 1593 | if(t->GetEMCALcluster()<=0 || t->GetEMCALcluster()>nclus) |
| | 1594 | continue; |
| | 1595 | AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(t->GetEMCALcluster())); |
| | 1596 | if(!c) |
| | 1597 | continue; |
| | 1598 | fEoverPvsE->Fill(c->E(),c->E()/t->P()); |
| | 1599 | } |
| | 1600 | for(int ic=0;ic<nclus;ic++){ |
| | 1601 | AliVCluster *c = dynamic_cast<AliVCluster*>(clusters->At(ic)); |
| | 1602 | //AliESDCaloCluster *c = (AliESDCaloCluster*)clusters->At(ic); |
| | 1603 | if(!c) |
| | 1604 | continue; |
| | 1605 | if(!c->IsEMCAL()) |
| | 1606 | continue; |
| | 1607 | Float_t clsPos[3] = {0,0,0}; |
| | 1608 | c->GetPosition(clsPos); |
| | 1609 | TVector3 clsVec(clsPos); |
| | 1610 | clsVec -= fVecPv; |
| | 1611 | Double_t cphi = clsVec.Phi(); |
| | 1612 | Double_t ceta = clsVec.Eta(); |
| | 1613 | Short_t id = -1; |
| | 1614 | GetMaxCellEnergy( c, id); |
| | 1615 | fEmcClusEClusCuts->Fill(c->E(),0); |
| | 1616 | fEmcClusEPhi->Fill(c->E(), cphi); |
| | 1617 | fEmcClusEEta->Fill(c->E(), ceta); |
| | 1618 | if(fMaxEClus>fECut){ |
| | 1619 | fEmcClusEPhiCut->Fill(c->E(), cphi); |
| | 1620 | fEmcClusEEtaCut->Fill(c->E(), ceta); |
| | 1621 | } |
| | 1622 | Double_t maxt=0; |
| | 1623 | if(fESDCells) |
| | 1624 | maxt = fESDCells->GetCellTime(id); |
| | 1625 | else if(fAODCells) |
| | 1626 | maxt = fAODCells->GetCellTime(id); |
| | 1627 | if(IsExotic(c)) |
| | 1628 | continue; |
| | 1629 | fEmcClusNotExo->Fill(c->E()); |
| | 1630 | fEmcClusEClusCuts->Fill(c->E(),1); |
| | 1631 | if(fClusIdFromTracks.Contains(Form("%d",ic))){ |
| | 1632 | fEmcClusETM2->Fill(c->E()); |
| | 1633 | fDetaDphiFromTM->Fill(c->GetTrackDz(),c->GetTrackDx()); |
| | 1634 | } |
| | 1635 | if(TMath::Abs(c->GetTrackDx())<0.03 && TMath::Abs(c->GetTrackDz())<0.02){ |
| | 1636 | fEmcClusETM1->Fill(c->E()); |
| | 1637 | continue; |
| | 1638 | } |
| | 1639 | fEmcClusEClusCuts->Fill(c->E(),2); |
| | 1640 | if(TMath::Abs(maxt)>30e-9) |
| | 1641 | continue; |
| | 1642 | fEmcClusEClusCuts->Fill(c->E(),3); |
| | 1643 | if(c->GetM02()>0.1) |
| | 1644 | fEmcClusEClusCuts->Fill(c->E(),4); |
| | 1645 | } |
| | 1646 | } |
| | 1647 | //________________________________________________________________________ |
| | 1648 | Double_t AliAnalysisTaskEMCALIsoPhoton::GetTrackMatchedPt(Int_t matchIndex) |
| | 1649 | { |
| | 1650 | Double_t pt = 0; |
| | 1651 | if(!fTracks) |
| | 1652 | return pt; |
| | 1653 | if(matchIndex<0 || matchIndex>fTracks->GetEntries()){ |
| | 1654 | if(fDebug) |
| | 1655 | printf("track-matched index out of track array range!!!\n"); |
| | 1656 | return pt; |
| | 1657 | } |
| | 1658 | AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(matchIndex)); |
| | 1659 | if(!track){ |
| | 1660 | if(fDebug) |
| | 1661 | printf("track-matched pointer does not exist!!!\n"); |
| | 1662 | return pt; |
| | 1663 | } |
| | 1664 | if(fESD){ |
| | 1665 | if(fPrTrCuts && fPrTrCuts->IsSelected(track)) |
| | 1666 | pt = track->Pt(); |
| | 1667 | else { |
| | 1668 | if(fCompTrCuts && fCompTrCuts->IsSelected(track)) |
| | 1669 | pt = track->Pt(); |
| | 1670 | else |
| | 1671 | return pt; |
| | 1672 | } |
| | 1673 | } |
| | 1674 | return pt; |
| | 1675 | } |
| | 1676 | //________________________________________________________________________ |
| | 1677 | void AliAnalysisTaskEMCALIsoPhoton::LoopOnCells() |
| | 1678 | { |
| | 1679 | AliVCaloCells *cells = 0; |
| | 1680 | cells = fESDCells; |
| | 1681 | if (!cells) |
| | 1682 | cells = fAODCells; |
| | 1683 | if(!cells) |
| | 1684 | return; |
| | 1685 | Double_t maxe = 0; |
| | 1686 | Double_t maxphi = -10; |
| | 1687 | Int_t ncells = cells->GetNumberOfCells(); |
| | 1688 | Double_t eta,phi; |
| | 1689 | for (Int_t i=0; i<ncells; i++) { |
| | 1690 | Short_t absid = TMath::Abs(cells->GetCellNumber(i)); |
| | 1691 | Double_t e = cells->GetCellAmplitude(absid); |
| | 1692 | if(e>0.05) |
| | 1693 | fNCells50++; |
| | 1694 | else |
| | 1695 | continue; |
| | 1696 | fGeom->EtaPhiFromIndex(absid,eta,phi); |
| | 1697 | if(maxe<e){ |
| | 1698 | maxe = e; |
| | 1699 | maxphi = phi; |
| | 1700 | } |
| | 1701 | } |
| | 1702 | fMaxCellEPhi->Fill(maxe,maxphi); |
| | 1703 | |
| | 1704 | } |
| | 1705 | //________________________________________________________________________ |
| | 1706 | bool AliAnalysisTaskEMCALIsoPhoton::IsExotic(AliVCluster *c) |
| | 1707 | { |
| | 1708 | bool isExo = 0; |
| | 1709 | Short_t id = -1; |
| | 1710 | Double_t Emax = GetMaxCellEnergy( c, id); |
| | 1711 | Double_t Ecross = GetCrossEnergy( c, id); |
| | 1712 | if((1-Ecross/Emax)>fExoticCut) |
| | 1713 | isExo = 1; |
| | 1714 | return isExo; |
| | 1715 | } |
| | 1716 | //________________________________________________________________________ |
| | 1717 | void AliAnalysisTaskEMCALIsoPhoton::Terminate(Option_t *) |
| | 1718 | { |
| | 1719 | // Called once at the end of the query. |
| | 1720 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental