| 3a66195e |
1 | #include "AliAnalysisTaskEMCALMesonGGSDM.h" |
| 2 | |
| 3 | #include <vector> |
| 4 | #include <Riostream.h> |
| 5 | #include <TChain.h> |
| 6 | #include <TTree.h> |
| 7 | #include <TF1.h> |
| 8 | #include <TH1F.h> |
| 9 | #include <TH2F.h> |
| 10 | #include <TH3F.h> |
| 11 | #include <TH1D.h> |
| 12 | #include <TH2D.h> |
| 13 | #include <TH3D.h> |
| 14 | #include <TCanvas.h> |
| 15 | #include <TList.h> |
| 16 | #include <TFile.h> |
| 17 | #include <TLorentzVector.h> |
| 18 | #include <TNtuple.h> |
| 19 | #include <TRandom3.h> |
| 20 | |
| 21 | #include "AliAnalysisTaskSE.h" |
| 22 | #include "AliAnalysisManager.h" |
| 23 | #include "AliStack.h" |
| 24 | #include "AliESDtrackCuts.h" |
| 25 | #include "AliESDEvent.h" |
| 26 | #include "AliESDInputHandler.h" |
| 27 | #include "AliAODEvent.h" |
| 28 | #include "AliMCEvent.h" |
| 29 | #include "AliEMCALGeometry.h" |
| 30 | #include "AliInputEventHandler.h" |
| 31 | #include "AliESDInputHandler.h" |
| 32 | #include "AliAODInputHandler.h" |
| 33 | |
| 34 | #include "AliEMCALRecoUtils.h" |
| 35 | #include "AliExternalTrackParam.h" |
| 36 | |
| 37 | // ROOT includes |
| 38 | #include <TGeoManager.h> |
| 39 | #include <TGeoMatrix.h> |
| 40 | #include <TGeoBBox.h> |
| 41 | #include <TH2F.h> |
| 42 | #include <TArrayI.h> |
| 43 | #include <TArrayF.h> |
| 44 | #include <TObjArray.h> |
| 45 | |
| 46 | // STEER includes |
| 47 | #include "AliVCluster.h" |
| 48 | #include "AliVCaloCells.h" |
| 49 | #include "AliLog.h" |
| 50 | #include "AliPID.h" |
| 51 | #include "AliESDEvent.h" |
| 52 | #include "AliAODEvent.h" |
| 53 | #include "AliESDtrack.h" |
| 54 | #include "AliAODTrack.h" |
| 55 | #include "AliExternalTrackParam.h" |
| 56 | #include "AliESDfriendTrack.h" |
| 57 | #include "AliTrackerBase.h" |
| 58 | |
| 59 | // EMCAL includes |
| 60 | #include "AliEMCALRecoUtils.h" |
| 61 | #include "AliEMCALGeometry.h" |
| 62 | #include "AliTrackerBase.h" |
| 63 | #include "AliEMCALCalibTimeDepCorrection.h" // Run dependent |
| 64 | #include "AliEMCALPIDUtils.h" |
| 65 | |
| 66 | #include "AliGenCocktailEventHeader.h" |
| 67 | |
| 2942f542 |
68 | using std::cout; |
| 69 | using std::endl; |
| 70 | |
| 3a66195e |
71 | ClassImp(AliAnalysisTaskEMCALMesonGGSDM) |
| 72 | |
| 73 | //________________________________________________________________________ |
| 74 | AliAnalysisTaskEMCALMesonGGSDM::AliAnalysisTaskEMCALMesonGGSDM() : |
| 75 | AliAnalysisTaskSE(), |
| 76 | fOutput(0), |
| 77 | fMcMode(0), |
| 78 | fMyMCType(0), |
| 79 | fRecalibrator(0), |
| 80 | fdRmin_ClustTrack(0), |
| 81 | fPhimin(0), |
| 82 | fPhimax(0), |
| 83 | fEtamin(0), |
| 84 | fEtamax(0), |
| 85 | fTrackCuts(0), |
| 86 | fEsdEv(0), |
| 87 | fAodEv(0), |
| 88 | h1_nClusters(0), |
| 89 | h1_zvtx(0), |
| 90 | h1_trigger(0), |
| 91 | h1_M(0), |
| 92 | h1_M_mix(0), |
| 93 | h1_E(0), |
| 94 | h2_PhiEtaCluster(0), |
| 95 | h2_PhiEtaClusterCut(0), |
| 96 | h2_PhiEtaMaxCell(0), |
| 97 | h2_PhiEtaMaxCellCut(0), |
| 98 | h1_dR_ClustTrk(0), |
| 99 | h2_gE_RecTruth(0), |
| 100 | h2_eop_E(0), |
| 101 | h2_eop_pT(0), |
| 102 | h2_E_time(0), |
| 103 | h1_Pi0TruthPt(0), |
| 104 | h1_K0Pi0TruthPt(0), |
| 105 | h1_PriPi0TruthPt(0), |
| 106 | h1_PhysPi0TruthPt(0), |
| 107 | h1_Pi0TruthPtEmcal(0), |
| 108 | h1_K0Pi0TruthPtEmcal(0), |
| 109 | h1_PriPi0TruthPtEmcal(0), |
| 110 | h1_PhysPi0TruthPtEmcal(0), |
| 111 | h1_Pi0TruthPtPhi2piEta065(0), |
| 112 | h1_K0Pi0TruthPtPhi2piEta065(0), |
| 113 | h1_PriPi0TruthPtPhi2piEta065(0), |
| 114 | h1_PhysPi0TruthPtPhi2piEta065(0), |
| 115 | h1_Pi0TruthPtPhi2piEta1(0), |
| 116 | h1_K0Pi0TruthPtPhi2piEta1(0), |
| 117 | h1_PriPi0TruthPtPhi2piEta1(0), |
| 118 | h1_PhysPi0TruthPtPhi2piEta1(0), |
| 119 | h2_Pi0TruthPhiEta(0), |
| 120 | h2_PriPi0TruthPhiEta(0), |
| 121 | h2_Pi0TruthPhiEtaEmcal(0), |
| 122 | h2_PriPi0TruthPhiEtaEmcal(0), |
| 123 | h1_TruthPhotonsEmcal(0), |
| 124 | h2_TruthPhotonsPhiEta(0), |
| 125 | h1_PhotonsEmcal(0), |
| 126 | h1_PhotonsNCellsCut(0), |
| 127 | h1_PhotonsTrackMatchCut(0), |
| 128 | h1_PhotonsAllCut(0), |
| 129 | h2_PhotonsPhiEtaIsEmcal(0), |
| 130 | h1_dR_RealMC(0), |
| 131 | h2_Mpt_Pri(0), |
| 132 | h2_Mpt_Sec(0), |
| 133 | h3_MptR_Sec(0), |
| 134 | h3_MptR_K0s(0), |
| 135 | h3_MptR_Mat(0), |
| 136 | h2_PtR_MatM(0), |
| 137 | h2_Mpt_Pri_conv(0), |
| 138 | h2_Mpt_Sec_conv(0), |
| 139 | h3_MptR_Sec_conv(0), |
| 140 | h3_MptR_K0s_conv(0), |
| 141 | h3_MptR_Mat_conv(0), |
| 142 | h1_eConversionR(0), |
| 143 | h1_PriPi0Mother(0), |
| 144 | h1_SecPi0Mother(0), |
| 145 | h1_Chi2(0), |
| 146 | h1_nTrkMatch(0), |
| 147 | h1_nCells(0), |
| 148 | h1_ClusterDisp(0), |
| 149 | h2_Ellipse(0), |
| 150 | h2_EtaPt(0), |
| 151 | //h2_Mpt(0), |
| 152 | h3_MptAsymm(0), |
| 153 | //h2_Mpt_mix(0), |
| 154 | h3_MptAsymm_mix(0), |
| 155 | h2_dphi_deta(0), |
| 156 | h2_dphi_deta_mix(0), |
| 157 | h2_DispRes(0), |
| 158 | h2_cells_M02(0), |
| 159 | TriggerList(0), |
| 160 | fHelperClass(0) |
| 161 | { |
| 162 | // Dummy constructor ALWAYS needed for I/O. |
| 163 | } |
| 164 | |
| 165 | //________________________________________________________________________ |
| 166 | AliAnalysisTaskEMCALMesonGGSDM::AliAnalysisTaskEMCALMesonGGSDM(const char *name) : |
| 167 | AliAnalysisTaskSE(name), |
| 168 | fOutput(0), |
| 169 | fMcMode(0), |
| 170 | fMyMCType(0), |
| 171 | fRecalibrator(0), |
| 172 | fdRmin_ClustTrack(0), |
| 173 | fPhimin(0), |
| 174 | fPhimax(0), |
| 175 | fEtamin(0), |
| 176 | fEtamax(0), |
| 177 | fTrackCuts(0), |
| 178 | fEsdEv(0), |
| 179 | fAodEv(0), |
| 180 | h1_nClusters(0), |
| 181 | h1_zvtx(0), |
| 182 | h1_trigger(0), |
| 183 | h1_M(0), |
| 184 | h1_M_mix(0), |
| 185 | h1_E(0), |
| 186 | h2_PhiEtaCluster(0), |
| 187 | h2_PhiEtaClusterCut(0), |
| 188 | h2_PhiEtaMaxCell(0), |
| 189 | h2_PhiEtaMaxCellCut(0), |
| 190 | h1_dR_ClustTrk(0), |
| 191 | h2_gE_RecTruth(0), |
| 192 | h2_eop_E(0), |
| 193 | h2_eop_pT(0), |
| 194 | h2_E_time(0), |
| 195 | h1_Pi0TruthPt(0), |
| 196 | h1_K0Pi0TruthPt(0), |
| 197 | h1_PriPi0TruthPt(0), |
| 198 | h1_PhysPi0TruthPt(0), |
| 199 | h1_Pi0TruthPtEmcal(0), |
| 200 | h1_K0Pi0TruthPtEmcal(0), |
| 201 | h1_PriPi0TruthPtEmcal(0), |
| 202 | h1_PhysPi0TruthPtEmcal(0), |
| 203 | h1_Pi0TruthPtPhi2piEta065(0), |
| 204 | h1_K0Pi0TruthPtPhi2piEta065(0), |
| 205 | h1_PriPi0TruthPtPhi2piEta065(0), |
| 206 | h1_PhysPi0TruthPtPhi2piEta065(0), |
| 207 | h1_Pi0TruthPtPhi2piEta1(0), |
| 208 | h1_K0Pi0TruthPtPhi2piEta1(0), |
| 209 | h1_PriPi0TruthPtPhi2piEta1(0), |
| 210 | h1_PhysPi0TruthPtPhi2piEta1(0), |
| 211 | h2_Pi0TruthPhiEta(0), |
| 212 | h2_PriPi0TruthPhiEta(0), |
| 213 | h2_Pi0TruthPhiEtaEmcal(0), |
| 214 | h2_PriPi0TruthPhiEtaEmcal(0), |
| 215 | h1_TruthPhotonsEmcal(0), |
| 216 | h2_TruthPhotonsPhiEta(0), |
| 217 | h1_PhotonsEmcal(0), |
| 218 | h1_PhotonsNCellsCut(0), |
| 219 | h1_PhotonsTrackMatchCut(0), |
| 220 | h1_PhotonsAllCut(0), |
| 221 | h2_PhotonsPhiEtaIsEmcal(0), |
| 222 | h1_dR_RealMC(0), |
| 223 | h2_Mpt_Pri(0), |
| 224 | h2_Mpt_Sec(0), |
| 225 | h3_MptR_Sec(0), |
| 226 | h3_MptR_K0s(0), |
| 227 | h3_MptR_Mat(0), |
| 228 | h2_PtR_MatM(0), |
| 229 | h2_Mpt_Pri_conv(0), |
| 230 | h2_Mpt_Sec_conv(0), |
| 231 | h3_MptR_Sec_conv(0), |
| 232 | h3_MptR_K0s_conv(0), |
| 233 | h3_MptR_Mat_conv(0), |
| 234 | h1_eConversionR(0), |
| 235 | h1_PriPi0Mother(0), |
| 236 | h1_SecPi0Mother(0), |
| 237 | h1_Chi2(0), |
| 238 | h1_nTrkMatch(0), |
| 239 | h1_nCells(0), |
| 240 | h1_ClusterDisp(0), |
| 241 | h2_Ellipse(0), |
| 242 | h2_EtaPt(0), |
| 243 | //h2_Mpt(0), |
| 244 | h3_MptAsymm(0), |
| 245 | //h2_Mpt_mix(0), |
| 246 | h3_MptAsymm_mix(0), |
| 247 | h2_dphi_deta(0), |
| 248 | h2_dphi_deta_mix(0), |
| 249 | h2_DispRes(0), |
| 250 | h2_cells_M02(0), |
| 251 | TriggerList(0), |
| 252 | fHelperClass(0) |
| 253 | { |
| 254 | // Constructor |
| 255 | // Define input and output slots here (never in the dummy constructor) |
| 256 | // Input slot #0 works with a TChain - it is connected to the default input container |
| 257 | // Output slot #1 writes into a TH1 container |
| 258 | |
| 259 | DefineOutput(1, TList::Class()); // for output list |
| 260 | } |
| 261 | |
| 262 | //________________________________________________________________________ |
| 263 | AliAnalysisTaskEMCALMesonGGSDM::~AliAnalysisTaskEMCALMesonGGSDM() |
| 264 | { |
| 265 | // Destructor. Clean-up the output list, but not the histograms that are put inside |
| 266 | // (the list is owner and will clean-up these histograms). Protect in PROOF case. |
| 267 | if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { |
| 268 | delete fOutput; |
| 269 | } |
| 270 | delete fTrackCuts; |
| 271 | } |
| 272 | |
| 273 | //________________________________________________________________________ |
| 274 | void AliAnalysisTaskEMCALMesonGGSDM::UserCreateOutputObjects() |
| 275 | { |
| 276 | // Create histograms |
| 277 | // Called once (on the worker node) |
| 278 | |
| 279 | fOutput = new TList(); |
| 280 | fOutput->SetOwner(); // IMPORTANT! |
| 281 | |
| 282 | fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); |
| 283 | |
| 284 | cout << "__________AliAnalysisTaskEMCALMesonGGSDM: Input settings__________" << endl; |
| 285 | cout << " fMcMode: " << fMcMode << endl; |
| 286 | cout << " fMyMCType: " << fMyMCType << endl; |
| 287 | cout << " fRecalibrator: " << fRecalibrator << endl; |
| 288 | cout << " dRmin_ClustTrack: " << fdRmin_ClustTrack << endl; |
| 289 | cout << " phi range: " << fPhimin << ", " << fPhimax << endl; |
| 290 | cout << " eta range: " << fEtamin << ", " << fEtamax << endl; |
| 291 | cout << " number of zvtx bins: " << zvtx_bins << endl; |
| 292 | cout << " number of mult bins: " << mult_bins << endl; |
| 293 | cout << " poolDepth: " << poolDepth << endl; |
| 294 | cout << endl; |
| 295 | |
| 296 | |
| 297 | //AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager(); |
| 298 | //AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler()); |
| 299 | //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse(); |
| 300 | // |
| 301 | //fPIDCombined = new AliPIDCombined(); |
| 302 | //fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES); |
| 303 | //fPIDCombined->SetDetectorMask(AliPIDResponse::kDetEMCAL); |
| 304 | //fPIDCombined->SetEnablePriors(kFALSE); |
| 305 | |
| 306 | double TotalNBins = 0.0; |
| 307 | |
| 308 | // Create histograms |
| 309 | Int_t nClustersbins = 501; |
| 310 | Float_t nClusterslow = -0.5, nClustersup = 500.5; |
| 311 | h1_nClusters = new TH1F("h1_nClusters", "# of clusters", nClustersbins, nClusterslow, nClustersup); |
| 312 | h1_nClusters->GetXaxis()->SetTitle("number of clusters/evt"); |
| 313 | h1_nClusters->GetYaxis()->SetTitle("counts"); |
| 314 | h1_nClusters->SetMarkerStyle(kFullCircle); |
| 315 | TotalNBins+=nClustersbins; |
| 316 | |
| 317 | Int_t nZvertexbins = 501; |
| 318 | Float_t Zvertexlow = -50.0, Zvertexup = 50.0; |
| 319 | h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup); |
| 320 | h1_zvtx->GetXaxis()->SetTitle("z_{vertex}"); |
| 321 | h1_zvtx->GetYaxis()->SetTitle("counts"); |
| 322 | h1_zvtx->SetMarkerStyle(kFullCircle); |
| 323 | TotalNBins+=nZvertexbins; |
| 324 | |
| 325 | h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5); |
| 326 | TotalNBins+=1001; |
| 327 | |
| 328 | Int_t Mbins = 3000; |
| 329 | Float_t Mlow = 0.0, Mup = 3.0; |
| 330 | h1_M = new TH1F("h1_M", "Invariant Mass", Mbins, Mlow, Mup); |
| 331 | h1_M->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 332 | h1_M->GetYaxis()->SetTitle("counts"); |
| 333 | h1_M->SetMarkerStyle(kFullCircle); |
| 334 | TotalNBins+=Mbins; |
| 335 | |
| 336 | h1_M_mix = new TH1F("h1_M_mix", "Invariant Mass (mixed events)", Mbins, Mlow, Mup); |
| 337 | h1_M_mix->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 338 | h1_M_mix->GetYaxis()->SetTitle("counts"); |
| 339 | h1_M_mix->SetMarkerStyle(kFullCircle); |
| 340 | TotalNBins+=Mbins; |
| 341 | |
| 342 | Int_t ptbins = 2000; |
| 343 | Float_t ptlow = 0.0, ptup = 20.0; |
| 344 | Int_t Ebins = 1000; |
| 345 | Float_t Elow = 0.0, Eup = 20.0; |
| 346 | h1_E = new TH1F("h1_E", "Cluster Energy in EMCal", Ebins, Elow, Eup); |
| 347 | h1_E->GetXaxis()->SetTitle("E [GeV]"); |
| 348 | h1_E->GetYaxis()->SetTitle("counts"); |
| 349 | h1_E->SetMarkerStyle(kFullCircle); |
| 350 | TotalNBins+=Ebins; |
| 351 | |
| 352 | h2_PhiEtaCluster = new TH2F("h2_PhiEtaCluster", "cluster phi vs eta", 400,1.362,3.178, 300,-0.728,0.728); |
| 353 | h2_PhiEtaCluster->GetXaxis()->SetTitle("#phi [rad]"); |
| 354 | h2_PhiEtaCluster->GetYaxis()->SetTitle("#eta"); |
| 355 | h2_PhiEtaCluster->GetZaxis()->SetTitle("hits"); |
| 356 | h2_PhiEtaCluster->SetMarkerStyle(kFullCircle); |
| 357 | TotalNBins+=400*300; |
| 358 | |
| 359 | h2_PhiEtaClusterCut = new TH2F("h2_PhiEtaClusterCut", "cluster phi vs eta (after cuts)", 400,1.362,3.178, 300,-0.728,0.728); |
| 360 | h2_PhiEtaClusterCut->GetXaxis()->SetTitle("#phi [rad]"); |
| 361 | h2_PhiEtaClusterCut->GetYaxis()->SetTitle("#eta"); |
| 362 | h2_PhiEtaClusterCut->GetZaxis()->SetTitle("hits"); |
| 363 | h2_PhiEtaClusterCut->SetMarkerStyle(kFullCircle); |
| 364 | TotalNBins+=400*300; |
| 365 | |
| 366 | // eta binning |
| 367 | Double_t EtaBins[97] = {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465}; |
| 368 | |
| 369 | // phi binning |
| 370 | Double_t PhiBins[125] = {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732, 1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128}; |
| 371 | |
| 372 | h2_PhiEtaMaxCell = new TH2F("h2_PhiEtaMaxCell", "maxcell phi vs eta", 124,PhiBins, 96,EtaBins); |
| 373 | h2_PhiEtaMaxCell->GetXaxis()->SetTitle("#phi [rad]"); |
| 374 | h2_PhiEtaMaxCell->GetYaxis()->SetTitle("#eta"); |
| 375 | h2_PhiEtaMaxCell->GetZaxis()->SetTitle("hits"); |
| 376 | h2_PhiEtaMaxCell->SetMarkerStyle(kFullCircle); |
| 377 | TotalNBins+=96*124; |
| 378 | |
| 379 | h2_PhiEtaMaxCellCut = new TH2F("h2_PhiEtaMaxCellCut", "maxcell phi vs eta (after cuts)", 124,PhiBins, 96,EtaBins); |
| 380 | h2_PhiEtaMaxCellCut->GetXaxis()->SetTitle("#phi [rad]"); |
| 381 | h2_PhiEtaMaxCellCut->GetYaxis()->SetTitle("#eta"); |
| 382 | h2_PhiEtaMaxCellCut->GetZaxis()->SetTitle("hits"); |
| 383 | h2_PhiEtaMaxCellCut->SetMarkerStyle(kFullCircle); |
| 384 | TotalNBins+=96*124; |
| 385 | |
| 386 | h1_dR_ClustTrk = new TH1F("h1_dR_ClustTrk", "Cluster-Track matching", 5000, -0.01, 5); |
| 387 | h1_dR_ClustTrk->GetXaxis()->SetTitle("dR [sqrt(d#phi^{2}+d#eta^{2})]"); |
| 388 | h1_dR_ClustTrk->GetYaxis()->SetTitle("N"); |
| 389 | h1_dR_ClustTrk->SetMarkerStyle(kFullCircle); |
| 390 | TotalNBins+=5000; |
| 391 | |
| 392 | h2_gE_RecTruth = new TH2F("h2_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 500,0,2); |
| 393 | h2_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]"); |
| 394 | h2_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}"); |
| 395 | h2_gE_RecTruth->GetZaxis()->SetTitle("counts"); |
| 396 | h2_gE_RecTruth->SetMarkerStyle(kFullCircle); |
| 397 | TotalNBins+=Ebins*500; |
| 398 | |
| 399 | h2_eop_E = new TH2F("h2_eop_E","E/p vs E (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3); |
| 400 | h2_eop_E->GetXaxis()->SetTitle("cluster Energy [GeV]"); |
| 401 | h2_eop_E->GetYaxis()->SetTitle("E/p"); |
| 402 | TotalNBins+=Ebins*1200; |
| 403 | |
| 404 | h2_eop_pT = new TH2F("h2_eop_pT","E/p vs p_{T} (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3); |
| 405 | h2_eop_pT->GetXaxis()->SetTitle("cluster Energy [GeV]"); |
| 406 | h2_eop_pT->GetYaxis()->SetTitle("E/p"); |
| 407 | TotalNBins+=Ebins*1200; |
| 408 | |
| 409 | h2_E_time = new TH2F("h2_E_time","cluster energy vs time", Ebins, Elow, Eup, 1000,-1e-6,1e-6); |
| 410 | h2_E_time->GetXaxis()->SetTitle("cluster Energy [GeV]"); |
| 411 | h2_E_time->GetYaxis()->SetTitle("time [s]"); |
| 412 | TotalNBins+=Ebins*1000; |
| 413 | |
| 414 | h1_Pi0TruthPt = new TH1F("h1_Pi0TruthPt", "P_{T} distribution for Truth Pi0's", ptbins, ptlow, ptup); |
| 415 | h1_Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 416 | h1_Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 417 | h1_Pi0TruthPt->SetMarkerStyle(kFullCircle); |
| 418 | TotalNBins+=ptbins; |
| 419 | |
| 420 | h1_K0Pi0TruthPt = new TH1F("h1_K0Pi0TruthPt", "P_{T} distribution for Truth Pi0's from K^{0}_{s} decays", ptbins, ptlow, ptup); |
| 421 | h1_K0Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 422 | h1_K0Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 423 | h1_K0Pi0TruthPt->SetMarkerStyle(kFullCircle); |
| 424 | TotalNBins+=ptbins; |
| 425 | |
| 426 | h1_PriPi0TruthPt = new TH1F("h1_PriPi0TruthPt", "P_{T} distribution for Truth Primary Pi0's", ptbins, ptlow, ptup); |
| 427 | h1_PriPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 428 | h1_PriPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 429 | h1_PriPi0TruthPt->SetMarkerStyle(kFullCircle); |
| 430 | TotalNBins+=ptbins; |
| 431 | |
| 432 | h1_PhysPi0TruthPt = new TH1F("h1_PhysPi0TruthPt", "P_{T} distribution for Truth Physical Primary Pi0's", ptbins, ptlow, ptup); |
| 433 | h1_PhysPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 434 | h1_PhysPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 435 | h1_PhysPi0TruthPt->SetMarkerStyle(kFullCircle); |
| 436 | TotalNBins+=ptbins; |
| 437 | |
| 438 | h1_Pi0TruthPtEmcal = new TH1F("h1_Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's (hit EMCal)", ptbins, ptlow, ptup); |
| 439 | h1_Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 440 | h1_Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 441 | h1_Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle); |
| 442 | TotalNBins+=ptbins; |
| 443 | |
| 444 | h1_K0Pi0TruthPtEmcal = new TH1F("h1_K0Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's from K^{0}_{s} decays (hit EMCal)", ptbins, ptlow, ptup); |
| 445 | h1_K0Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 446 | h1_K0Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 447 | h1_K0Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle); |
| 448 | TotalNBins+=ptbins; |
| 449 | |
| 450 | h1_PriPi0TruthPtEmcal = new TH1F("h1_PriPi0TruthPtEmcal", "P_{T} distribution for Truth Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup); |
| 451 | h1_PriPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 452 | h1_PriPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 453 | h1_PriPi0TruthPtEmcal->SetMarkerStyle(kFullCircle); |
| 454 | TotalNBins+=ptbins; |
| 455 | |
| 456 | h1_PhysPi0TruthPtEmcal = new TH1F("h1_PhysPi0TruthPtEmcal", "P_{T} distribution for Truth Physical Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup); |
| 457 | h1_PhysPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 458 | h1_PhysPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 459 | h1_PhysPi0TruthPtEmcal->SetMarkerStyle(kFullCircle); |
| 460 | TotalNBins+=ptbins; |
| 461 | |
| 462 | h1_Pi0TruthPtPhi2piEta065 = new TH1F("h1_Pi0TruthPtPhi2piEta065", |
| 463 | "P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 464 | h1_Pi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 465 | h1_Pi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 466 | h1_Pi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle); |
| 467 | TotalNBins+=ptbins; |
| 468 | |
| 469 | h1_K0Pi0TruthPtPhi2piEta065 = new TH1F("h1_K0Pi0TruthPtPhi2piEta065", |
| 470 | "P_{T} for Truth Pi0's (K^{0}_{s} decays) [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 471 | h1_K0Pi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 472 | h1_K0Pi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 473 | h1_K0Pi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle); |
| 474 | TotalNBins+=ptbins; |
| 475 | |
| 476 | h1_PriPi0TruthPtPhi2piEta065 = new TH1F("h1_PriPi0TruthPtPhi2piEta065", |
| 477 | "P_{T} for Primary Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 478 | h1_PriPi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 479 | h1_PriPi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 480 | h1_PriPi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle); |
| 481 | TotalNBins+=ptbins; |
| 482 | |
| 483 | h1_PhysPi0TruthPtPhi2piEta065 = new TH1F("h1_PhysPi0TruthPtPhi2piEta065", |
| 484 | "P_{T} for Truth Pi0's (not from material) [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 485 | h1_PhysPi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 486 | h1_PhysPi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 487 | h1_PhysPi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle); |
| 488 | TotalNBins+=ptbins; |
| 489 | |
| 490 | h1_Pi0TruthPtPhi2piEta1 = new TH1F("h1_Pi0TruthPtPhi2piEta1", |
| 491 | "P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 492 | h1_Pi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 493 | h1_Pi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 494 | h1_Pi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle); |
| 495 | TotalNBins+=ptbins; |
| 496 | |
| 497 | h1_K0Pi0TruthPtPhi2piEta1 = new TH1F("h1_K0Pi0TruthPtPhi2piEta1", |
| 498 | "P_{T} for Truth Pi0's (k^{0}_{s} decays) [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 499 | h1_K0Pi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 500 | h1_K0Pi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 501 | h1_K0Pi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle); |
| 502 | TotalNBins+=ptbins; |
| 503 | |
| 504 | h1_PriPi0TruthPtPhi2piEta1 = new TH1F("h1_PriPi0TruthPtPhi2piEta1", |
| 505 | "P_{T} for Primary Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 506 | h1_PriPi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 507 | h1_PriPi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 508 | h1_PriPi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle); |
| 509 | TotalNBins+=ptbins; |
| 510 | |
| 511 | h1_PhysPi0TruthPtPhi2piEta1 = new TH1F("h1_PhysPi0TruthPtPhi2piEta1", |
| 512 | "P_{T} for Truth Pi0's (not from material) [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup); |
| 513 | h1_PhysPi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 514 | h1_PhysPi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 515 | h1_PhysPi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle); |
| 516 | TotalNBins+=ptbins; |
| 517 | |
| 518 | h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10); |
| 519 | h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]"); |
| 520 | h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta "); |
| 521 | TotalNBins+=380*200; |
| 522 | |
| 523 | h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10); |
| 524 | h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]"); |
| 525 | h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta "); |
| 526 | TotalNBins+=380*200; |
| 527 | |
| 528 | h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5); |
| 529 | h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]"); |
| 530 | h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta "); |
| 531 | TotalNBins+=380*150; |
| 532 | |
| 533 | h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5); |
| 534 | h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]"); |
| 535 | h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta "); |
| 536 | TotalNBins+=380*150; |
| 537 | |
| 538 | h1_TruthPhotonsEmcal = new TH1F("h1_TruthPhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup); |
| 539 | h1_TruthPhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 540 | h1_TruthPhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 541 | h1_TruthPhotonsEmcal->SetMarkerStyle(kFullCircle); |
| 542 | TotalNBins+=ptbins; |
| 543 | |
| 544 | h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta", |
| 545 | "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5); |
| 546 | h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]"); |
| 547 | h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta "); |
| 548 | h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle); |
| 549 | TotalNBins+=380*150; |
| 550 | |
| 551 | h1_PhotonsEmcal = new TH1F("h1_PhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup); |
| 552 | h1_PhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 553 | h1_PhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 554 | h1_PhotonsEmcal->SetMarkerStyle(kFullCircle); |
| 555 | TotalNBins+=ptbins; |
| 556 | |
| 557 | h1_PhotonsNCellsCut = new TH1F("h1_PhotonsNCellsCut", "P_{T} distribution for #gamma's that survive NCells cut", ptbins, ptlow, ptup); |
| 558 | h1_PhotonsNCellsCut->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 559 | h1_PhotonsNCellsCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 560 | h1_PhotonsNCellsCut->SetMarkerStyle(kFullCircle); |
| 561 | TotalNBins+=ptbins; |
| 562 | |
| 563 | h1_PhotonsTrackMatchCut = new TH1F("h1_PhotonsTrackMatchCut", "P_{T} distribution for #gamma's that survive TrackMatch cut", ptbins, ptlow, ptup); |
| 564 | h1_PhotonsTrackMatchCut->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 565 | h1_PhotonsTrackMatchCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 566 | h1_PhotonsTrackMatchCut->SetMarkerStyle(kFullCircle); |
| 567 | TotalNBins+=ptbins; |
| 568 | |
| 569 | h1_PhotonsAllCut = new TH1F("h1_PhotonsAllCut", "P_{T} distribution for #gamma's that survive All cut", ptbins, ptlow, ptup); |
| 570 | h1_PhotonsAllCut->GetXaxis()->SetTitle("P_{T} (GeV/c)"); |
| 571 | h1_PhotonsAllCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)"); |
| 572 | h1_PhotonsAllCut->SetMarkerStyle(kFullCircle); |
| 573 | TotalNBins+=ptbins; |
| 574 | |
| 575 | h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal", |
| 576 | "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 100,-1.0,1.0); |
| 577 | h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]"); |
| 578 | h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta "); |
| 579 | TotalNBins+=380*100; |
| 580 | |
| 581 | h1_dR_RealMC = new TH1F("h1_dR_RealMC", "P_{T} distribution for #gamma's that survive All cut", 2000, -0.01, 10); |
| 582 | h1_dR_RealMC->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})"); |
| 583 | h1_dR_RealMC->GetYaxis()->SetTitle("N"); |
| 584 | h1_dR_RealMC->SetMarkerStyle(kFullCircle); |
| 585 | TotalNBins+=2000; |
| 586 | |
| 587 | h2_Mpt_Pri = new TH2F("h2_Mpt_Pri", "mass vs pT for primary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup); |
| 588 | h2_Mpt_Pri->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 589 | h2_Mpt_Pri->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 590 | h2_Mpt_Pri->SetMarkerStyle(kFullCircle); |
| 591 | TotalNBins+=Mbins*ptbins; |
| 592 | |
| 593 | h2_Mpt_Sec = new TH2F("h2_Mpt_Sec", "mass vs pT for secondary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup); |
| 594 | h2_Mpt_Sec->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 595 | h2_Mpt_Sec->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 596 | h2_Mpt_Sec->SetMarkerStyle(kFullCircle); |
| 597 | TotalNBins+=Mbins*ptbins; |
| 598 | |
| 599 | h3_MptR_Sec = new TH3F("h3_MptR_Sec", "mass vs pT vs production radius for secondary pions", 500,0,0.5, 100,0,20, 300,0,600); |
| 600 | h3_MptR_Sec->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 601 | h3_MptR_Sec->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 602 | h3_MptR_Sec->GetZaxis()->SetTitle("production radius [cm]"); |
| 603 | h3_MptR_Sec->SetMarkerStyle(kFullCircle); |
| 604 | TotalNBins+=500*100*300; |
| 605 | |
| 606 | h3_MptR_K0s = new TH3F("h3_MptR_K0s", "mass vs pT vs production radius for K0s pions", 500,0,0.5, 100,0,20, 300,0,600); |
| 607 | h3_MptR_K0s->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 608 | h3_MptR_K0s->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 609 | h3_MptR_K0s->GetZaxis()->SetTitle("production radius [cm]"); |
| 610 | h3_MptR_K0s->SetMarkerStyle(kFullCircle); |
| 611 | TotalNBins+=500*100*300; |
| 612 | |
| 613 | h3_MptR_Mat = new TH3F("h3_MptR_Mat", "mass vs pT vs production radius for material pions", 500,0,0.5, 100,0,20, 300,0,600); |
| 614 | h3_MptR_Mat->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 615 | h3_MptR_Mat->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 616 | h3_MptR_Mat->GetZaxis()->SetTitle("production radius [cm]"); |
| 617 | h3_MptR_Mat->SetMarkerStyle(kFullCircle); |
| 618 | TotalNBins+=500*100*300; |
| 619 | |
| 620 | h2_PtR_MatM = new TH2F("h2_PtR_MatM", "pT vs production radius for merged material pions (pi mass assumed)", 100,0,20, 300,0,600); |
| 621 | h2_PtR_MatM->GetXaxis()->SetTitle("p_{T} [GeV/c]"); |
| 622 | h2_PtR_MatM->GetYaxis()->SetTitle("production radius [cm]"); |
| 623 | h2_PtR_MatM->SetMarkerStyle(kFullCircle); |
| 624 | TotalNBins+=100*300; |
| 625 | |
| 626 | h2_Mpt_Pri_conv = new TH2F("h2_Mpt_Pri_conv", "mass vs pT for primary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup); |
| 627 | h2_Mpt_Pri_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 628 | h2_Mpt_Pri_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 629 | h2_Mpt_Pri_conv->SetMarkerStyle(kFullCircle); |
| 630 | TotalNBins+=Mbins*ptbins; |
| 631 | |
| 632 | h2_Mpt_Sec_conv = new TH2F("h2_Mpt_Sec_conv", "mass vs pT for secondary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup); |
| 633 | h2_Mpt_Sec_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 634 | h2_Mpt_Sec_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 635 | h2_Mpt_Sec_conv->SetMarkerStyle(kFullCircle); |
| 636 | TotalNBins+=Mbins*ptbins; |
| 637 | |
| 638 | h3_MptR_Sec_conv = new TH3F("h3_MptR_Sec_conv", "mass vs pT vs production radius for secondary pions", 500,0,0.5, 100,0,20, 300,0,600); |
| 639 | h3_MptR_Sec_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 640 | h3_MptR_Sec_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 641 | h3_MptR_Sec_conv->GetZaxis()->SetTitle("production radius [cm]"); |
| 642 | h3_MptR_Sec_conv->SetMarkerStyle(kFullCircle); |
| 643 | TotalNBins+=500*100*300; |
| 644 | |
| 645 | h3_MptR_K0s_conv = new TH3F("h3_MptR_K0s_conv", "mass vs pT vs production radius for K0s pions", 500,0,0.5, 100,0,20, 300,0,600); |
| 646 | h3_MptR_K0s_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 647 | h3_MptR_K0s_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 648 | h3_MptR_K0s_conv->GetZaxis()->SetTitle("production radius [cm]"); |
| 649 | h3_MptR_K0s_conv->SetMarkerStyle(kFullCircle); |
| 650 | TotalNBins+=500*100*300; |
| 651 | |
| 652 | h3_MptR_Mat_conv = new TH3F("h3_MptR_Mat_conv", "mass vs pT vs production radius for material pions", 500,0,0.5, 100,0,20, 300,0,600); |
| 653 | h3_MptR_Mat_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 654 | h3_MptR_Mat_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 655 | h3_MptR_Mat_conv->GetZaxis()->SetTitle("production radius [cm]"); |
| 656 | h3_MptR_Mat_conv->SetMarkerStyle(kFullCircle); |
| 657 | TotalNBins+=500*100*300; |
| 658 | |
| 659 | h1_eConversionR = new TH1F("h1_eConversionR", "conversion point (radius)", 600,0,600); |
| 660 | h1_eConversionR->GetXaxis()->SetTitle("production radius [cm]"); |
| 661 | h1_eConversionR->SetMarkerStyle(kFullCircle); |
| 662 | TotalNBins+=600; |
| 663 | |
| 664 | h1_PriPi0Mother = new TH1F("h1_PriPi0Mother", "primary pi0 mother ID", 12001,-6000.5,6000.5); |
| 665 | h1_PriPi0Mother->GetXaxis()->SetTitle("#pi^{0} mother ID"); |
| 666 | h1_PriPi0Mother->SetMarkerStyle(kFullCircle); |
| 667 | TotalNBins+=12001; |
| 668 | |
| 669 | h1_SecPi0Mother = new TH1F("h1_SecPi0Mother", "secondray pi0 mother ID", 12001,-6000.5,6000.5); |
| 670 | h1_SecPi0Mother->GetXaxis()->SetTitle("#pi^{0} mother ID"); |
| 671 | h1_SecPi0Mother->SetMarkerStyle(kFullCircle); |
| 672 | TotalNBins+=12001; |
| 673 | |
| 674 | Int_t chi2bins = 100; |
| 675 | Float_t chi2low = -2, chi2up = 2; |
| 676 | h1_Chi2 = new TH1F("h1_Chi2","#chi^{2} distribution for reconstructed",chi2bins, chi2low, chi2up); |
| 677 | h1_Chi2->GetXaxis()->SetTitle("#chi^{2}"); |
| 678 | h1_Chi2->GetYaxis()->SetTitle("counts"); |
| 679 | TotalNBins+=chi2bins; |
| 680 | |
| 681 | h1_nTrkMatch = new TH1F("h1_nTrkMatch","number of matched tracks",14, -1.5, 5.5); |
| 682 | h1_nTrkMatch->GetXaxis()->SetTitle("nTracksMatched"); |
| 683 | h1_nTrkMatch->GetYaxis()->SetTitle("counts"); |
| 684 | TotalNBins+=14; |
| 685 | |
| 686 | h1_ClusterDisp = new TH1F("h1_ClusterDisp","Dispersion of CaloCluster",1000, -1, 3); |
| 687 | h1_ClusterDisp->GetXaxis()->SetTitle("cluster->GetClusterDisp()"); |
| 688 | h1_ClusterDisp->GetYaxis()->SetTitle("counts"); |
| 689 | TotalNBins+=1000; |
| 690 | |
| 691 | h2_Ellipse = new TH2F("h2_Ellipse","Ellipse axis M20 vs M02",500, -0.01, 1, 500, -0.01, 1); |
| 692 | h2_Ellipse->GetXaxis()->SetTitle("cluster->GetM20()"); |
| 693 | h2_Ellipse->GetYaxis()->SetTitle("cluster->GetM02()"); |
| 694 | h2_Ellipse->GetZaxis()->SetTitle("counts"); |
| 695 | TotalNBins+=500*500; |
| 696 | |
| 697 | Int_t etabins = 150; |
| 698 | Float_t etalow = -1.5, etaup = 1.5; |
| 699 | h2_EtaPt = new TH2F("h2_EtaPt","Cluster Energy vs ",etabins, etalow, etaup, ptbins, ptlow, ptup); |
| 700 | h2_EtaPt->GetXaxis()->SetTitle("E [GeV]"); |
| 701 | h2_EtaPt->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 702 | TotalNBins+=etabins*ptbins; |
| 703 | |
| 704 | h3_MptAsymm = new TH3F("h3_MptAsymm","mass vs p_{T} vs Asymm cut",Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5); |
| 705 | h3_MptAsymm->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 706 | h3_MptAsymm->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 707 | h3_MptAsymm->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)"); |
| 708 | TotalNBins+=Mbins*ptbins*3.0; |
| 709 | |
| 710 | h3_MptAsymm_mix = new TH3F("h3_MptAsymm_mix","mass vs p_{T} vs Asymm cut (mixed events)",Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5); |
| 711 | h3_MptAsymm_mix->GetXaxis()->SetTitle("mass [GeV/c^{2}]"); |
| 712 | h3_MptAsymm_mix->GetYaxis()->SetTitle("p_{T} [GeV/c]"); |
| 713 | h3_MptAsymm_mix->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)"); |
| 714 | TotalNBins+=Mbins*ptbins*3.0; |
| 715 | |
| 716 | h2_dphi_deta = new TH2F("h2_dphi_deta","#Delta#phi vs #Delta#eta", 349,-1.5,5, 400,-2.0,2.0); |
| 717 | h2_dphi_deta->GetXaxis()->SetTitle("#Delta#phi"); |
| 718 | h2_dphi_deta->GetYaxis()->SetTitle("#Delta#eta"); |
| 719 | TotalNBins+=349*400; |
| 720 | |
| 721 | h2_dphi_deta_mix = new TH2F("h2_dphi_deta_mix","#Delta#phi vs #Delta#eta (mixed events)", 349,-1.5,5, 400,-2.0,2.0); |
| 722 | h2_dphi_deta_mix->GetXaxis()->SetTitle("#Delta#phi"); |
| 723 | h2_dphi_deta_mix->GetYaxis()->SetTitle("#Delta#eta"); |
| 724 | TotalNBins+=349*400; |
| 725 | |
| 726 | h2_DispRes = new TH2F("h2_DispRes", "zvtx info", 500,-0.01,1, 500,-0.1,2); |
| 727 | h2_DispRes->GetXaxis()->SetTitle("EvtVtx->GetDispersion()"); |
| 728 | h2_DispRes->GetYaxis()->SetTitle("EvtVtx->GetZRes()"); |
| 729 | h2_DispRes->GetZaxis()->SetTitle("counts"); |
| 730 | TotalNBins+=500*500; |
| 731 | |
| 732 | h2_cells_M02 = new TH2F("h2_cells_M02", "nCells vs M02", 204,-1.5,100.5, 500,-1,1.5); |
| 733 | h2_cells_M02->GetXaxis()->SetTitle("nCells"); |
| 734 | h2_cells_M02->GetYaxis()->SetTitle("M02"); |
| 735 | h2_cells_M02->GetZaxis()->SetTitle("counts"); |
| 736 | TotalNBins+=204*500; |
| 737 | |
| 738 | cout << endl << "Total number of bins in booked histograms: " << TotalNBins << endl << endl; |
| 739 | |
| 740 | // Initialize helper class (for vertex selection & pile up correction) |
| 741 | fHelperClass = new AliAnalysisUtils(); |
| 742 | |
| 743 | //TFile *f = OpenFile(1); |
| 744 | //TDirectory::TContext context(f); |
| 745 | |
| 746 | fOutput->Add(h1_nClusters); |
| 747 | fOutput->Add(h1_zvtx); |
| 748 | fOutput->Add(h1_trigger); |
| 749 | fOutput->Add(h1_M); |
| 750 | fOutput->Add(h1_M_mix); |
| 751 | fOutput->Add(h1_E); |
| 752 | fOutput->Add(h2_PhiEtaCluster); |
| 753 | fOutput->Add(h2_PhiEtaClusterCut); |
| 754 | fOutput->Add(h2_PhiEtaMaxCell); |
| 755 | fOutput->Add(h2_PhiEtaMaxCellCut); |
| 756 | fOutput->Add(h1_dR_ClustTrk); |
| 757 | fOutput->Add(h2_gE_RecTruth); |
| 758 | fOutput->Add(h2_eop_E); |
| 759 | fOutput->Add(h2_eop_pT); |
| 760 | fOutput->Add(h2_E_time); |
| 761 | fOutput->Add(h1_Pi0TruthPt); |
| 762 | fOutput->Add(h1_K0Pi0TruthPt); |
| 763 | fOutput->Add(h1_PriPi0TruthPt); |
| 764 | fOutput->Add(h1_PhysPi0TruthPt); |
| 765 | fOutput->Add(h1_Pi0TruthPtEmcal); |
| 766 | fOutput->Add(h1_K0Pi0TruthPtEmcal); |
| 767 | fOutput->Add(h1_PriPi0TruthPtEmcal); |
| 768 | fOutput->Add(h1_PhysPi0TruthPtEmcal); |
| 769 | fOutput->Add(h1_Pi0TruthPtPhi2piEta065); |
| 770 | fOutput->Add(h1_K0Pi0TruthPtPhi2piEta065); |
| 771 | fOutput->Add(h1_PriPi0TruthPtPhi2piEta065); |
| 772 | fOutput->Add(h1_PhysPi0TruthPtPhi2piEta065); |
| 773 | fOutput->Add(h1_Pi0TruthPtPhi2piEta1); |
| 774 | fOutput->Add(h1_K0Pi0TruthPtPhi2piEta1); |
| 775 | fOutput->Add(h1_PriPi0TruthPtPhi2piEta1); |
| 776 | fOutput->Add(h1_PhysPi0TruthPtPhi2piEta1); |
| 777 | fOutput->Add(h2_Pi0TruthPhiEta); |
| 778 | fOutput->Add(h2_PriPi0TruthPhiEta); |
| 779 | fOutput->Add(h2_Pi0TruthPhiEtaEmcal); |
| 780 | fOutput->Add(h2_PriPi0TruthPhiEtaEmcal); |
| 781 | fOutput->Add(h1_TruthPhotonsEmcal); |
| 782 | fOutput->Add(h2_TruthPhotonsPhiEta); |
| 783 | fOutput->Add(h1_PhotonsEmcal); |
| 784 | fOutput->Add(h1_PhotonsNCellsCut); |
| 785 | fOutput->Add(h1_PhotonsTrackMatchCut); |
| 786 | fOutput->Add(h1_PhotonsAllCut); |
| 787 | fOutput->Add(h2_PhotonsPhiEtaIsEmcal); |
| 788 | fOutput->Add(h1_dR_RealMC); |
| 789 | fOutput->Add(h2_Mpt_Pri); |
| 790 | fOutput->Add(h2_Mpt_Sec); |
| 791 | fOutput->Add(h3_MptR_Sec); |
| 792 | fOutput->Add(h3_MptR_K0s); |
| 793 | fOutput->Add(h3_MptR_Mat); |
| 794 | fOutput->Add(h2_PtR_MatM); |
| 795 | fOutput->Add(h2_Mpt_Pri_conv); |
| 796 | fOutput->Add(h2_Mpt_Sec_conv); |
| 797 | fOutput->Add(h3_MptR_Sec_conv); |
| 798 | fOutput->Add(h3_MptR_K0s_conv); |
| 799 | fOutput->Add(h3_MptR_Mat_conv); |
| 800 | fOutput->Add(h1_eConversionR); |
| 801 | fOutput->Add(h1_PriPi0Mother); |
| 802 | fOutput->Add(h1_SecPi0Mother); |
| 803 | fOutput->Add(h1_Chi2); |
| 804 | fOutput->Add(h1_nTrkMatch); |
| 805 | fOutput->Add(h1_ClusterDisp); |
| 806 | fOutput->Add(h2_Ellipse); |
| 807 | fOutput->Add(h2_EtaPt); |
| 808 | fOutput->Add(h3_MptAsymm); |
| 809 | fOutput->Add(h3_MptAsymm_mix); |
| 810 | fOutput->Add(h2_dphi_deta); |
| 811 | fOutput->Add(h2_dphi_deta_mix); |
| 812 | fOutput->Add(h2_DispRes); |
| 813 | fOutput->Add(h2_cells_M02); |
| 814 | |
| 815 | // Post data for ALL output slots >0 here, |
| 816 | // To get at least an empty histogram |
| 817 | // 1 is the outputnumber of a certain weg of task 1 |
| 818 | PostData(1, fOutput); |
| 819 | } |
| 820 | |
| 821 | //________________________________________________________________________ |
| 822 | void AliAnalysisTaskEMCALMesonGGSDM::UserExec(Option_t *) |
| 823 | { |
| 824 | // Main loop Called for each event |
| 825 | |
| 826 | AliMCEvent *mcEvent = MCEvent(); |
| 827 | Bool_t isMC = bool(mcEvent);//is this the right way to do this? |
| 828 | |
| 829 | TRandom3 randy; randy.SetSeed(0); |
| 830 | unsigned int iskip = -1; |
| 831 | TLorentzVector ParentMix; |
| 832 | |
| 833 | double recalScale = 1.0; |
| 834 | |
| 835 | AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager(); |
| 836 | |
| 837 | AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler()); |
| 838 | AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler()); |
| 839 | if (!aodH && !esdH) Printf("ERROR: Could not get ESD or AODInputHandler"); |
| 840 | |
| 841 | if(esdH) fEsdEv = esdH->GetEvent(); |
| 842 | else if(aodH) fAodEv = aodH->GetEvent(); |
| 843 | else{ |
| 844 | AliFatal("Neither ESD nor AOD event found"); |
| 845 | return; |
| 846 | } |
| 847 | |
| 848 | // get pointer to reconstructed event |
| 849 | AliVEvent *event = InputEvent(); |
| 850 | if (!event){ |
| 851 | AliError("Pointer == 0, this can not happen!"); return;} |
| 852 | //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event); |
| 853 | //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event); |
| 854 | //if (!fEsdEv){ |
| 855 | //AliError("Cannot get the ESD event"); return;} |
| 856 | |
| 857 | /* |
| 858 | fHelperClass->SetCutOnZVertexSPD(kFALSE);//does the zvtx have to match the spd vertex? |
| 859 | fHelperClass->SetMaxVtxZ(1.0e6);//i set this myself later.. |
| 860 | // simply makes sure that there is at least 1 contributer to the zvtx determination. |
| 861 | // this should only remove the *extra* events at zvtx==0. |
| 862 | if(!fHelperClass->IsVertexSelected2013pA(event)) |
| 863 | return; |
| 864 | */ |
| 865 | |
| 866 | Int_t iTrigger = 0; |
| 867 | if (fEsdEv) iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs(); |
| 868 | else if (fAodEv) iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs(); |
| 869 | |
| 870 | char saythis[500]; |
| 871 | Int_t iTriggerBin = 0; |
| 872 | for(unsigned long j=0; j<TriggerList.size(); j++){ |
| 873 | if(iTrigger==TriggerList[j]) |
| 874 | iTriggerBin=j+1; |
| 875 | } |
| 876 | if(iTriggerBin==0){ |
| 877 | TriggerList.push_back(iTrigger); |
| 878 | iTriggerBin=TriggerList.size(); |
| 879 | } |
| 880 | |
| 881 | h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1); |
| 882 | sprintf(saythis,"%d",iTrigger); |
| 883 | h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis); |
| 884 | |
| 885 | if(fEsdEv){ |
| 886 | TString trigClasses = fEsdEv->GetFiredTriggerClasses(); |
| 887 | // remove "fast cluster events": |
| 888 | if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL")) |
| 889 | return; |
| 890 | } |
| 891 | else if(fAodEv){ |
| 892 | TString trigClasses = fAodEv->GetFiredTriggerClasses(); |
| 893 | // remove "fast cluster events": |
| 894 | if (trigClasses.Contains("FAST") && !trigClasses.Contains("ALL")) |
| 895 | return; |
| 896 | } |
| 897 | |
| 898 | if (fEsdEv){ |
| 899 | if(!(fEsdEv->GetPrimaryVertex()->GetStatus())) return; |
| 900 | } |
| 901 | //else if (fAodEv){ |
| 902 | //if(!(fAodEv->GetPrimaryVertex()->GetStatus())) return; |
| 903 | //} |
| 904 | |
| 905 | Double_t vertDisp=0.0; |
| 906 | Double_t vertZres=0.0; |
| 907 | Bool_t vertIsfromZ=0; |
| 908 | if (fEsdEv){ |
| 909 | vertDisp = fEsdEv->GetPrimaryVertex()->GetDispersion(); |
| 910 | vertZres = fEsdEv->GetPrimaryVertex()->GetZRes(); |
| 911 | vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ(); |
| 912 | } |
| 913 | else if (fAodEv){ |
| 914 | vertDisp = 0; |
| 915 | vertZres = 0; |
| 916 | vertIsfromZ = 0; |
| 917 | } |
| 918 | |
| 919 | h2_DispRes->Fill(vertDisp, vertZres); |
| 920 | // if vertex is from spd vertexZ, require more stringent cut |
| 921 | if (vertIsfromZ) { |
| 922 | if (vertDisp>0.02 || vertZres>0.25 ) |
| 923 | return; // bad vertex from VertexerZ |
| 924 | } |
| 925 | |
| 926 | |
| 927 | Int_t nclusters=0; |
| 928 | if(fEsdEv){ |
| 929 | //Int_t evtN = fEsdEv->GetEventNumberInFile(); |
| 930 | //Int_t ntracks = fEsdEv->GetNumberOfTracks(); |
| 931 | nclusters = fEsdEv->GetNumberOfCaloClusters(); |
| 932 | } |
| 933 | else if(fAodEv){ |
| 934 | //Int_t evtN = fAodEv->GetEventNumberInFile(); |
| 935 | //Int_t ntracks = fAodEv->GetNumberOfTracks(); |
| 936 | nclusters = fAodEv->GetNumberOfCaloClusters(); |
| 937 | } |
| 938 | |
| 939 | // EMCal cluster loop for reconstructed event |
| 940 | //numberofclusters set above! |
| 941 | TLorentzVector Photon1, Photon2, Parent; |
| 942 | Double_t vertex[3]; |
| 943 | Double_t E1=0.0; |
| 944 | Double_t E2=0.0; |
| 945 | Double_t vertZ=0.0; |
| 946 | if (fEsdEv) vertZ = fEsdEv->GetPrimaryVertex()->GetZ(); |
| 947 | else if (fAodEv) vertZ = fAodEv->GetPrimaryVertex()->GetZ(); |
| 948 | |
| 949 | h1_zvtx->Fill(vertZ); |
| 950 | //zvertex cut: |
| 951 | if(fabs(vertZ)>10.0) |
| 952 | return; |
| 953 | |
| 954 | h1_nClusters->Fill(nclusters); |
| 955 | |
| 956 | int izvtx = GetZvtxBin(vertZ); |
| 957 | int imult = GetMultBin(nclusters); |
| 958 | |
| 959 | //cout << iskip << " " << izvtx << " " << imult << endl; |
| 960 | //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl; |
| 961 | |
| 962 | |
| 963 | |
| 964 | //######################### ~~~~~~~~~~~ ################################## |
| 965 | //######################### STARTING MC ################################## |
| 966 | //######################### ~~~~~~~~~~~ ################################## |
| 967 | |
| 968 | if(isMC){ |
| 969 | int isPrimary = 0; |
| 970 | int isMaterialSec = 0; |
| 971 | int isK0sDecay = 0; |
| 972 | |
| 973 | if (!mcEvent){ |
| 974 | cout << "no MC event" << endl; |
| 975 | return; |
| 976 | } |
| 977 | |
| 978 | const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex(); |
| 979 | if (!evtVtx) |
| 980 | return; |
| 981 | |
| 982 | mcEvent->PreReadAll(); |
| 983 | |
| 984 | Int_t nTracksMC = mcEvent->GetNumberOfTracks(); |
| 985 | Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries(); |
| 986 | Int_t MyEMCPion = -1; |
| 987 | for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) { |
| 988 | AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack)); |
| 989 | if (!mcP) |
| 990 | continue; |
| 991 | |
| 992 | if(iTrack<nPTracksMC) isPrimary = 1; |
| 993 | else isPrimary = 0; |
| 994 | |
| 995 | isK0sDecay = 0; |
| 996 | if(mcP->GetMother()>-1){ |
| 997 | if( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 310 || |
| 998 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -310 ) |
| 999 | isK0sDecay = 1; |
| 1000 | } |
| 1001 | |
| 1002 | // it's a pion !! |
| 1003 | if(mcP->PdgCode() != 111) |
| 1004 | continue; |
| 1005 | |
| 1006 | MyEMCPion = 0; |
| 1007 | if(strcmp(fMyMCType,"ANY")==0) |
| 1008 | MyEMCPion = 1; |
| 1009 | else if(IsMyMCHeaderType(iTrack, fMyMCType, mcEvent)){ |
| 1010 | //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << " i: " << iTrack << " pdg: " << mcP->PdgCode() << " pT: " << mcP->Pt() << endl; |
| 1011 | //cout << "iTrack: " << iTrack << " nPrimaryMC: " << nPTracksMC << " nTracksMC: " << nTracksMC << endl; |
| 1012 | MyEMCPion = 1; |
| 1013 | } |
| 1014 | |
| 1015 | //if(MyEMCPion) |
| 1016 | //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << " i: " << iTrack << " pdg: " << mcP->PdgCode() << " pT: " << mcP->Pt() << endl; |
| 1017 | //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << " i: " << iTrack << " pdg: " << static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack+1))->PdgCode() << " pT: " << mcP->Pt() << endl; |
| 1018 | |
| 1019 | if(MyEMCPion!=1 && isPrimary==1) |
| 1020 | continue; |
| 1021 | |
| 1022 | |
| 1023 | if(isPrimary==1 && mcP->GetMother()>-1) |
| 1024 | h1_PriPi0Mother->Fill( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ); |
| 1025 | else if(isPrimary==0 && mcP->GetMother()>-1) |
| 1026 | h1_SecPi0Mother->Fill( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ); |
| 1027 | |
| 1028 | Int_t daughter[2] = {-1,-1}; |
| 1029 | daughter[0] = mcP->GetFirstDaughter(); |
| 1030 | daughter[1] = mcP->GetLastDaughter(); |
| 1031 | |
| 1032 | if (daughter[0]<0) continue; |
| 1033 | if (daughter[1]<0) daughter[1]=daughter[0]; |
| 1034 | if (daughter[1]-daughter[0] != 1) continue; |
| 1035 | |
| 1036 | Int_t eIndexofConvertedPhoton[2] = {-1,-1}; |
| 1037 | |
| 1038 | bool bacc = true; |
| 1039 | bool binp = true; |
| 1040 | Double_t eta_d[2] = {0.0,0.0}; |
| 1041 | Double_t phi_d[2] = {0.0,0.0}; |
| 1042 | for (Int_t daughter_index=0; daughter_index<2; daughter_index++){ |
| 1043 | const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(daughter[daughter_index])); |
| 1044 | eta_d[daughter_index] = dmc->Eta(); |
| 1045 | phi_d[daughter_index] = dmc->Phi(); |
| 1046 | if(!(dmc->PdgCode()==22)) binp = false; |
| 1047 | if(!(dmc->PdgCode()==22 && |
| 1048 | eta_d[daughter_index]>fEtamin && eta_d[daughter_index]<fEtamax && |
| 1049 | phi_d[daughter_index]>fPhimin && phi_d[daughter_index]<fPhimax)) bacc = false; |
| 1050 | |
| 1051 | if(dmc->GetFirstDaughter()>0 && dmc->GetLastDaughter()>0) { |
| 1052 | // get the photons's daughters... |
| 1053 | const AliMCParticle *dmcd1 = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dmc->GetFirstDaughter())); |
| 1054 | const AliMCParticle *dmcd2 = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dmc->GetLastDaughter())); |
| 1055 | Double_t productionR1 = TMath::Sqrt(dmcd1->Xv()*dmcd1->Xv() + dmcd1->Yv()*dmcd1->Yv()); |
| 1056 | if(bacc) h1_eConversionR->Fill(productionR1); |
| 1057 | // check if this is a conversion... |
| 1058 | if( (dmcd1->PdgCode()== -1.0*dmcd2->PdgCode()) && |
| 1059 | (dmcd1->PdgCode()==11 || dmcd1->PdgCode()==-11) && |
| 1060 | productionR1<440.0){ |
| 1061 | //find the conv e with highest energy, assign it to be that photon decay product. |
| 1062 | if( dmcd1->E() > dmcd2->E() ) |
| 1063 | eIndexofConvertedPhoton[daughter_index] = dmc->GetFirstDaughter(); |
| 1064 | else |
| 1065 | eIndexofConvertedPhoton[daughter_index] = dmc->GetLastDaughter(); |
| 1066 | } |
| 1067 | } |
| 1068 | } |
| 1069 | |
| 1070 | if(binp!=true) |
| 1071 | continue; |
| 1072 | |
| 1073 | // primary particle |
| 1074 | //Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) + |
| 1075 | // (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY())); |
| 1076 | //if(dR <= 0.01) isPrimary = 1; |
| 1077 | //else isPrimary = 0; |
| 1078 | |
| 1079 | isMaterialSec = 0; |
| 1080 | if(isPrimary!=1){ |
| 1081 | if(((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 2212 || //proton |
| 1082 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -2212 || //anti-proton |
| 1083 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 2112 || //neutron |
| 1084 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -2112 || //anti-neutron |
| 1085 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 321 || //K+ |
| 1086 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -321 || //K- |
| 1087 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == 211 || //pi+ |
| 1088 | ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -211 //pi- |
| 1089 | ) |
| 1090 | isMaterialSec = 1; |
| 1091 | } |
| 1092 | |
| 1093 | h1_Pi0TruthPt ->Fill(mcP->Pt()); |
| 1094 | if(isK0sDecay) h1_K0Pi0TruthPt->Fill(mcP->Pt()); |
| 1095 | h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta()); |
| 1096 | |
| 1097 | if(isPrimary==1){ |
| 1098 | h1_PriPi0TruthPt ->Fill(mcP->Pt()); |
| 1099 | h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta()); |
| 1100 | } |
| 1101 | if(isPrimary!=1 && isMaterialSec!=1) h1_PhysPi0TruthPt->Fill(mcP->Pt()); |
| 1102 | |
| 1103 | if(mcP->Eta()<-1.0 || mcP->Eta()>1.0) |
| 1104 | continue; |
| 1105 | |
| 1106 | h1_Pi0TruthPtPhi2piEta1 ->Fill(mcP->Pt()); |
| 1107 | if(isPrimary==1) |
| 1108 | h1_PriPi0TruthPtPhi2piEta1 ->Fill(mcP->Pt()); |
| 1109 | if(isK0sDecay) |
| 1110 | h1_K0Pi0TruthPtPhi2piEta1 ->Fill(mcP->Pt()); |
| 1111 | if(isPrimary!=1 && isMaterialSec!=1) |
| 1112 | h1_PhysPi0TruthPtPhi2piEta1 ->Fill(mcP->Pt()); |
| 1113 | |
| 1114 | if(mcP->Eta()>fEtamin && mcP->Eta()<fEtamax){ |
| 1115 | h1_Pi0TruthPtPhi2piEta065 ->Fill(mcP->Pt()); |
| 1116 | if(isPrimary==1) |
| 1117 | h1_PriPi0TruthPtPhi2piEta065 ->Fill(mcP->Pt()); |
| 1118 | if(isK0sDecay) |
| 1119 | h1_K0Pi0TruthPtPhi2piEta065 ->Fill(mcP->Pt()); |
| 1120 | if(isPrimary!=1 && isMaterialSec!=1) |
| 1121 | h1_PhysPi0TruthPtPhi2piEta065 ->Fill(mcP->Pt()); |
| 1122 | } |
| 1123 | |
| 1124 | |
| 1125 | if(binp && bacc){// 2 Photons hit the EMCAL! |
| 1126 | |
| 1127 | Int_t Nfoundphotons = 0; |
| 1128 | Int_t iFoundphotons[10] = {0,0,0,0,0,0,0,0,0,0}; |
| 1129 | Int_t Nfoundelectrons = 0; |
| 1130 | Int_t iFoundelectrons[10] = {0,0,0,0,0,0,0,0,0,0}; |
| 1131 | for (Int_t daughter_index=0; daughter_index<2; daughter_index++){//both truth photons. (also includes conversions..) |
| 1132 | |
| 1133 | const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(daughter[daughter_index])); |
| 1134 | |
| 1135 | h1_TruthPhotonsEmcal->Fill(dmc->Pt()); |
| 1136 | h2_TruthPhotonsPhiEta->Fill(dmc->Phi(),dmc->Eta()); |
| 1137 | |
| 1138 | for(int i=0; i<nclusters; i++) { |
| 1139 | Bool_t matches_pion_photon = 0; |
| 1140 | |
| 1141 | AliESDCaloCluster* esdCluster=NULL; |
| 1142 | AliAODCaloCluster* aodCluster=NULL; |
| 1143 | if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster |
| 1144 | else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster |
| 1145 | |
| 1146 | Double_t clustMC_phi, clustMC_eta; |
| 1147 | |
| 1148 | if(fEsdEv){ |
| 1149 | |
| 1150 | if(esdCluster->IsEMCAL()){ |
| 1151 | |
| 1152 | Float_t pos[3] = {0,0,0}; |
| 1153 | esdCluster->GetPosition(pos); |
| 1154 | TVector3 vpos(pos); |
| 1155 | //h1_Phi->Fill(vpos.Phi()); |
| 1156 | clustMC_phi = vpos.Phi(); |
| 1157 | clustMC_eta = vpos.Eta(); |
| 1158 | |
| 1159 | Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) + |
| 1160 | (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi)); |
| 1161 | h1_dR_RealMC->Fill(dR); |
| 1162 | matches_pion_photon = 0; |
| 1163 | //if(dR<=0.04) matches_pion_photon = 1; |
| 1164 | |
| 1165 | |
| 1166 | TArrayI *TruthLabelsA = esdCluster->GetLabelsArray(); |
| 1167 | if(TruthLabelsA){ |
| 1168 | Int_t trackindex = TruthLabelsA->At(0); |
| 1169 | if( trackindex==daughter[daughter_index] ){ |
| 1170 | matches_pion_photon = 1; |
| 1171 | iFoundphotons[Nfoundphotons] = i; |
| 1172 | Nfoundphotons++; |
| 1173 | } |
| 1174 | else if( trackindex==eIndexofConvertedPhoton[daughter_index] ){ |
| 1175 | iFoundelectrons[Nfoundelectrons] = i; |
| 1176 | Nfoundelectrons++; |
| 1177 | } |
| 1178 | AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(trackindex)); |
| 1179 | |
| 1180 | if(matches_pion_photon){ |
| 1181 | |
| 1182 | h1_PhotonsEmcal->Fill(esdCluster->E()); |
| 1183 | h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta); |
| 1184 | if(esdCluster->GetNCells()>=2) |
| 1185 | h1_PhotonsNCellsCut->Fill(esdCluster->E()); |
| 1186 | if(esdCluster->GetNTracksMatched()==0) |
| 1187 | h1_PhotonsTrackMatchCut->Fill(esdCluster->E()); |
| 1188 | if(esdCluster->GetNCells()>=2 && esdCluster->GetNTracksMatched()==0) |
| 1189 | h1_PhotonsAllCut->Fill(esdCluster->E()); |
| 1190 | |
| 1191 | if(esdCluster->GetNCells()>=2){ |
| 1192 | recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator); |
| 1193 | h2_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E())); |
| 1194 | } |
| 1195 | }//if(matches_pion_photon) |
| 1196 | |
| 1197 | }//if Truthlabels exists |
| 1198 | vpos.Delete(); |
| 1199 | }//if(IsEMCAL()) |
| 1200 | |
| 1201 | }//if(fEsdEv) |
| 1202 | else if(fAodEv){ |
| 1203 | |
| 1204 | if(aodCluster->IsEMCAL()){ |
| 1205 | |
| 1206 | Float_t pos[3] = {0,0,0}; |
| 1207 | aodCluster->GetPosition(pos); |
| 1208 | TVector3 vpos(pos); |
| 1209 | //h1_Phi->Fill(vpos.Phi()); |
| 1210 | clustMC_phi = vpos.Phi(); |
| 1211 | clustMC_eta = vpos.Eta(); |
| 1212 | |
| 1213 | Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) + |
| 1214 | (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi)); |
| 1215 | h1_dR_RealMC->Fill(dR); |
| 1216 | matches_pion_photon = 0; |
| 1217 | if(dR<=0.04) matches_pion_photon = 1; |
| 1218 | |
| 1219 | //TArrayI *TruthLabelsA = aodCluster->GetLabelsArray(); |
| 1220 | //if(TruthLabelsA){ |
| 1221 | // Int_t trackindex = TruthLabelsA->At(0); |
| 1222 | // if( trackindex==daughter[daughter_index] ) |
| 1223 | // matches_pion_photon = 1; |
| 1224 | // AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(trackindex)); |
| 1225 | |
| 1226 | if(matches_pion_photon){ |
| 1227 | |
| 1228 | h1_PhotonsEmcal->Fill(aodCluster->E()); |
| 1229 | h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta); |
| 1230 | if(aodCluster->GetNCells()>=2) |
| 1231 | h1_PhotonsNCellsCut->Fill(aodCluster->E()); |
| 1232 | if(aodCluster->GetNTracksMatched()==0) |
| 1233 | h1_PhotonsTrackMatchCut->Fill(aodCluster->E()); |
| 1234 | if(aodCluster->GetNCells()>=2 && aodCluster->GetNTracksMatched()==0) |
| 1235 | h1_PhotonsAllCut->Fill(aodCluster->E()); |
| 1236 | |
| 1237 | //if(aodCluster->GetNCells()>=2){ |
| 1238 | //recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator); |
| 1239 | //h2_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E())); |
| 1240 | //} |
| 1241 | }//if(matches_pion_photon) |
| 1242 | |
| 1243 | //}//if Truthlabels exists |
| 1244 | vpos.Delete(); |
| 1245 | }//if(IsEMCAL()) |
| 1246 | |
| 1247 | }//if(fAodEv) |
| 1248 | |
| 1249 | }//loop over nclusters. |
| 1250 | |
| 1251 | }//both truth photons. |
| 1252 | |
| 1253 | if(Nfoundphotons>1){ |
| 1254 | AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster |
| 1255 | AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundphotons[1]); // pointer to EMCal cluster |
| 1256 | |
| 1257 | if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){ |
| 1258 | |
| 1259 | recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator); |
| 1260 | E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ |
| 1261 | fEsdEv->GetVertex()->GetXYZ(vertex); |
| 1262 | esdCluster1->GetMomentum(Photon1,vertex); |
| 1263 | Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ |
| 1264 | Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ |
| 1265 | Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ |
| 1266 | |
| 1267 | recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator); |
| 1268 | E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ |
| 1269 | fEsdEv->GetVertex()->GetXYZ(vertex); |
| 1270 | esdCluster2->GetMomentum(Photon2,vertex); |
| 1271 | Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ |
| 1272 | Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ |
| 1273 | Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ |
| 1274 | |
| 1275 | Parent = TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2); |
| 1276 | |
| 1277 | //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() ); |
| 1278 | Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) + |
| 1279 | (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY())); |
| 1280 | if(isPrimary==1){ |
| 1281 | //cout << "Primary production vertex: " << productionR << endl; |
| 1282 | h2_Mpt_Pri->Fill(Parent.M(),Parent.Pt()); |
| 1283 | } |
| 1284 | else{ |
| 1285 | //cout << "Secondary production vertex: " << productionR << endl; |
| 1286 | h2_Mpt_Sec ->Fill(Parent.M(),Parent.Pt()); |
| 1287 | h3_MptR_Sec->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1288 | if(isK0sDecay) |
| 1289 | h3_MptR_K0s->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1290 | if(isMaterialSec) |
| 1291 | h3_MptR_Mat->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1292 | } |
| 1293 | }//both good clusters |
| 1294 | }//found 2 photons. |
| 1295 | else if(Nfoundphotons==1){ |
| 1296 | int mergedPion = 0; |
| 1297 | AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster |
| 1298 | |
| 1299 | TArrayI *TruthLabelsA = esdCluster1->GetLabelsArray(); |
| 1300 | if(TruthLabelsA){ |
| 1301 | if(TruthLabelsA->GetSize()>1){ |
| 1302 | Int_t trackindex[2]; |
| 1303 | trackindex[0] = TruthLabelsA->At(0); |
| 1304 | trackindex[1] = TruthLabelsA->At(1); |
| 1305 | if( (trackindex[0]==daughter[0] && trackindex[1]==daughter[1]) || |
| 1306 | (trackindex[0]==daughter[1] && trackindex[1]==daughter[0]) ){ |
| 1307 | mergedPion = 1; |
| 1308 | } |
| 1309 | if(mergedPion==1){ |
| 1310 | recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator); |
| 1311 | E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ |
| 1312 | Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) + |
| 1313 | (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY())); |
| 1314 | h2_PtR_MatM->Fill(E1,productionR); |
| 1315 | }//if merged pion. |
| 1316 | }//truthlabel.size > 1 |
| 1317 | }//if truthlabels |
| 1318 | }// Nfoundphotons==1 |
| 1319 | |
| 1320 | if(Nfoundphotons==1 && Nfoundelectrons==1){ |
| 1321 | AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster |
| 1322 | AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundelectrons[0]); // pointer to EMCal cluster |
| 1323 | |
| 1324 | if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){ |
| 1325 | |
| 1326 | recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator); |
| 1327 | E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ |
| 1328 | fEsdEv->GetVertex()->GetXYZ(vertex); |
| 1329 | esdCluster1->GetMomentum(Photon1,vertex); |
| 1330 | Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ |
| 1331 | Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ |
| 1332 | Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ |
| 1333 | |
| 1334 | recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator); |
| 1335 | E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ |
| 1336 | fEsdEv->GetVertex()->GetXYZ(vertex); |
| 1337 | esdCluster2->GetMomentum(Photon2,vertex); |
| 1338 | Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ |
| 1339 | Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ |
| 1340 | Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ |
| 1341 | |
| 1342 | Parent = TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2); |
| 1343 | |
| 1344 | //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() ); |
| 1345 | Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) + |
| 1346 | (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY())); |
| 1347 | if(isPrimary==1){ |
| 1348 | //cout << "Primary production vertex: " << productionR << endl; |
| 1349 | h2_Mpt_Pri_conv->Fill(Parent.M(),Parent.Pt()); |
| 1350 | } |
| 1351 | else{ |
| 1352 | //cout << "Secondary production vertex: " << productionR << endl; |
| 1353 | h2_Mpt_Sec_conv ->Fill(Parent.M(),Parent.Pt()); |
| 1354 | h3_MptR_Sec_conv->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1355 | if(isK0sDecay) |
| 1356 | h3_MptR_K0s_conv->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1357 | if(isMaterialSec) |
| 1358 | h3_MptR_Mat_conv->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1359 | } |
| 1360 | }//both good clusters |
| 1361 | }// Nfoundphotons==1 && Nfoundelectrons==1 |
| 1362 | else if(Nfoundelectrons==2){ |
| 1363 | AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundelectrons[0]); // pointer to EMCal cluster |
| 1364 | AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundelectrons[1]); // pointer to EMCal cluster |
| 1365 | |
| 1366 | if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){ |
| 1367 | |
| 1368 | recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator); |
| 1369 | E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ |
| 1370 | fEsdEv->GetVertex()->GetXYZ(vertex); |
| 1371 | esdCluster1->GetMomentum(Photon1,vertex); |
| 1372 | Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ |
| 1373 | Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ |
| 1374 | Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ |
| 1375 | |
| 1376 | recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator); |
| 1377 | E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ |
| 1378 | fEsdEv->GetVertex()->GetXYZ(vertex); |
| 1379 | esdCluster2->GetMomentum(Photon2,vertex); |
| 1380 | Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ |
| 1381 | Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ |
| 1382 | Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ |
| 1383 | |
| 1384 | Parent = TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2); |
| 1385 | |
| 1386 | //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() ); |
| 1387 | Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) + |
| 1388 | (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY())); |
| 1389 | if(isPrimary==1){ |
| 1390 | //cout << "Primary production vertex: " << productionR << endl; |
| 1391 | h2_Mpt_Pri_conv->Fill(Parent.M(),Parent.Pt()); |
| 1392 | } |
| 1393 | else{ |
| 1394 | //cout << "Secondary production vertex: " << productionR << endl; |
| 1395 | h2_Mpt_Sec_conv ->Fill(Parent.M(),Parent.Pt()); |
| 1396 | h3_MptR_Sec_conv->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1397 | if(isK0sDecay) |
| 1398 | h3_MptR_K0s_conv->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1399 | if(isMaterialSec) |
| 1400 | h3_MptR_Mat_conv->Fill(Parent.M(),Parent.Pt(),productionR); |
| 1401 | } |
| 1402 | }//both good clusters |
| 1403 | }// Nfoundelectrons==2 |
| 1404 | |
| 1405 | h1_Pi0TruthPtEmcal ->Fill(mcP->Pt()); |
| 1406 | if(isK0sDecay) |
| 1407 | h1_K0Pi0TruthPtEmcal ->Fill(mcP->Pt()); |
| 1408 | h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta()); |
| 1409 | |
| 1410 | if(isPrimary==1){ |
| 1411 | h1_PriPi0TruthPtEmcal ->Fill(mcP->Pt()); |
| 1412 | h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta()); |
| 1413 | } |
| 1414 | if(isPrimary!=1 && isMaterialSec!=1) h1_PhysPi0TruthPtEmcal->Fill(mcP->Pt()); |
| 1415 | |
| 1416 | |
| 1417 | }// 2 Photons hit the EMCAL! |
| 1418 | |
| 1419 | }//for(nTracksMC) |
| 1420 | |
| 1421 | }//if(isMC) |
| 1422 | |
| 1423 | //######################### ~~~~~~~~~~~~ ################################## |
| 1424 | //######################### DONE WITH MC ################################## |
| 1425 | //######################### ~~~~~~~~~~~~ ################################## |
| 1426 | |
| 1427 | |
| 1428 | for(int i=0; i<nclusters; i++) { |
| 1429 | |
| 1430 | AliESDCaloCluster* esdCluster=NULL; |
| 1431 | AliAODCaloCluster* aodCluster=NULL; |
| 1432 | if (fEsdEv) esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster |
| 1433 | else if (fAodEv) aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster |
| 1434 | if(!esdCluster && !aodCluster) { |
| 1435 | AliError(Form("ERROR: Could not retrieve any (ESD or AOD) Cluster %d",i)); |
| 1436 | continue; |
| 1437 | } |
| 1438 | |
| 1439 | if(fEsdEv){ |
| 1440 | |
| 1441 | recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator); |
| 1442 | |
| 1443 | //uncomment this to do the track matching (1 of 3 lines, esd part)!! |
| 1444 | //Bool_t MatchesToTrack = 0; |
| 1445 | if(esdCluster->IsEMCAL()){ |
| 1446 | |
| 1447 | Float_t pos[3] = {0,0,0}; |
| 1448 | Short_t maxCellID = -1; |
| 1449 | Float_t celleta, cellphi; |
| 1450 | esdCluster->GetPosition(pos); |
| 1451 | TVector3 clusterPosition(pos); |
| 1452 | h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta()); |
| 1453 | GetMaxCellEnergy(esdCluster, maxCellID); |
| 1454 | AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance(); |
| 1455 | fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi); |
| 1456 | h2_PhiEtaMaxCell->Fill(cellphi,celleta); |
| 1457 | |
| 1458 | // _______________Track loop for reconstructed event_____________ |
| 1459 | for(Int_t itrk = 0; itrk < fEsdEv->GetNumberOfTracks(); itrk++) { |
| 1460 | AliESDtrack* esdTrack = fEsdEv->GetTrack(itrk); // pointer to reconstructed to track |
| 1461 | if(!esdTrack) { |
| 1462 | AliError(Form("ERROR: Could not retrieve any (ESD) track %d",itrk)); |
| 1463 | continue; |
| 1464 | } |
| 1465 | |
| 1466 | Double_t posTrk[3] = {0,0,0}; |
| 1467 | esdTrack->GetXYZ(posTrk); |
| 1468 | TVector3 vposTrk(posTrk); |
| 1469 | |
| 1470 | Double_t fMass = 0.139; |
| 1471 | Double_t fStepSurface = 20.; |
| 1472 | Float_t etaproj, phiproj, pttrackproj; |
| 1473 | |
| 1474 | AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam()); |
| 1475 | if(!trackParam) continue; |
| 1476 | AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj); |
| 1477 | |
| 1478 | double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) + |
| 1479 | (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) ); |
| 1480 | |
| 1481 | h1_dR_ClustTrk->Fill(dR_clusttrk); |
| 1482 | |
| 1483 | //uncomment this to do the track matching (2 of 3 lines, esd part)!! |
| 1484 | //if(dR_clusttrk<fdRmin_ClustTrack) |
| 1485 | //MatchesToTrack = 1; |
| 1486 | |
| 1487 | }//_____________________________nTracks__________________________ |
| 1488 | |
| 1489 | h2_cells_M02 ->Fill(esdCluster->GetNCells(),esdCluster->GetM02()); |
| 1490 | h2_Ellipse ->Fill(esdCluster->GetM20(),esdCluster->GetM02()); |
| 1491 | h1_Chi2 ->Fill(esdCluster->Chi2());//always -1. |
| 1492 | h1_nTrkMatch ->Fill(esdCluster->GetNTracksMatched()); |
| 1493 | h1_ClusterDisp->Fill(esdCluster->GetDispersion()); |
| 1494 | h2_E_time ->Fill(esdCluster->E(),esdCluster->GetTOF()); |
| 1495 | |
| 1496 | TArrayI *TrackLabels = esdCluster->GetTracksMatched(); |
| 1497 | if(TrackLabels){ |
| 1498 | if(TrackLabels->GetSize()>0){ |
| 1499 | Int_t trackindex = TrackLabels->At(0); |
| 1500 | AliESDtrack* matchingT = fEsdEv->GetTrack(trackindex); // pointer to reconstructed to track |
| 1501 | |
| 1502 | recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator); |
| 1503 | h2_eop_E ->Fill(esdCluster->E()*recalScale, esdCluster->E()*recalScale/matchingT->P()); |
| 1504 | h2_eop_pT->Fill(matchingT->Pt(), esdCluster->E()*recalScale/matchingT->P()); |
| 1505 | } |
| 1506 | } |
| 1507 | |
| 1508 | //uncomment this to do the track matching (2 of 3 lines, esd part)!! |
| 1509 | //if(isGoodEsdCluster(esdCluster) && !MatchesToTrack){ |
| 1510 | if(isGoodEsdCluster(esdCluster)){ |
| 1511 | recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator); |
| 1512 | E1 = esdCluster->E()*recalScale;// TOTAL HACK - JJ |
| 1513 | fEsdEv->GetVertex()->GetXYZ(vertex); |
| 1514 | esdCluster->GetMomentum(Photon1,vertex); |
| 1515 | Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ |
| 1516 | Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ |
| 1517 | Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ |
| 1518 | Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) ); |
| 1519 | h1_E->Fill(E1); |
| 1520 | h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta()); |
| 1521 | h2_PhiEtaMaxCellCut->Fill(cellphi,celleta); |
| 1522 | } |
| 1523 | clusterPosition.Delete(); |
| 1524 | }//if(esdCluster->isEMCAL()) |
| 1525 | }//if(fEsdEv) |
| 1526 | else if(fAodEv){ |
| 1527 | |
| 1528 | recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator); |
| 1529 | |
| 1530 | //uncomment this to do the track matching (1 of 3 lines, aod part)!! |
| 1531 | //Bool_t MatchesToTrack = 0; |
| 1532 | if(aodCluster->IsEMCAL()){ |
| 1533 | |
| 1534 | Float_t pos[3] = {0,0,0}; |
| 1535 | Short_t maxCellID = -1; |
| 1536 | Float_t celleta, cellphi; |
| 1537 | aodCluster->GetPosition(pos); |
| 1538 | TVector3 clusterPosition(pos); |
| 1539 | h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta()); |
| 1540 | GetMaxCellEnergy(aodCluster, maxCellID); |
| 1541 | AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance(); |
| 1542 | fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi); |
| 1543 | h2_PhiEtaMaxCell->Fill(cellphi,celleta); |
| 1544 | |
| 1545 | // _______________Track loop for reconstructed event_____________ |
| 1546 | for(Int_t itrk = 0; itrk < fAodEv->GetNumberOfTracks(); itrk++) { |
| 1547 | AliAODTrack* aodTrack = fAodEv->GetTrack(itrk); // pointer to reconstructed to track |
| 1548 | if(!aodTrack) { |
| 1549 | AliError(Form("ERROR: Could not retrieve any (AOD) track %d",itrk)); |
| 1550 | continue; |
| 1551 | } |
| 1552 | |
| 1553 | Double_t posTrk[3] = {0,0,0}; |
| 1554 | Double_t momTrk[3] = {0,0,0}; |
| 1555 | aodTrack->GetXYZ(posTrk); |
| 1556 | aodTrack->GetPxPyPz(momTrk); |
| 1557 | //TVector3 vposTrk(posTrk); |
| 1558 | |
| 1559 | //#################################################################################################### |
| 1560 | // |
| 1561 | // commented all this stuff just to satisfy aliroot warnings. |
| 1562 | // but I may need it again if I want to do the track matching for aods. |
| 1563 | /* |
| 1564 | Double_t fMass = 0.139; |
| 1565 | Double_t fStepSurface = 20.; |
| 1566 | Float_t etaproj=0.0; |
| 1567 | Float_t phiproj=0.0; |
| 1568 | Float_t pttrackproj=0.0; |
| 1569 | |
| 1570 | Double_t cv[21] = {0.0}; |
| 1571 | aodTrack->GetCovarianceXYZPxPyPz(cv); |
| 1572 | AliExternalTrackParam *trackParam = new AliExternalTrackParam(posTrk,momTrk,cv,aodTrack->Charge()); |
| 1573 | //AliExternalTrackParam emcalParam(*trackParam); |
| 1574 | //AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(aodTrack->GetInnerParam()); |
| 1575 | if(!trackParam) continue; |
| 1576 | ////AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj); |
| 1577 | //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&emcalParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj); |
| 1578 | delete trackParam; |
| 1579 | |
| 1580 | //Constantin's implementation... gives funny result. |
| 1581 | //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(aodTrack,440.0); |
| 1582 | //phiproj = aodTrack->GetTrackPhiOnEMCal(); |
| 1583 | //etaproj = aodTrack->GetTrackPhiOnEMCal(); |
| 1584 | |
| 1585 | double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) + |
| 1586 | (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) ); |
| 1587 | |
| 1588 | h1_dR_ClustTrk->Fill(dR_clusttrk); |
| 1589 | */ |
| 1590 | //#################################################################################################### |
| 1591 | |
| 1592 | //uncomment this to do the track matching (2 of 3 lines, aod part)!! |
| 1593 | //if(dR_clusttrk<fdRmin_ClustTrack) |
| 1594 | //MatchesToTrack = 1; |
| 1595 | |
| 1596 | |
| 1597 | }//_____________________________nTracks__________________________ |
| 1598 | |
| 1599 | h2_cells_M02 ->Fill(aodCluster->GetNCells(),aodCluster->GetM02()); |
| 1600 | h2_Ellipse ->Fill(aodCluster->GetM20(),aodCluster->GetM02()); |
| 1601 | h1_Chi2 ->Fill(aodCluster->Chi2());//always -1. |
| 1602 | h1_nTrkMatch ->Fill(aodCluster->GetNTracksMatched()); |
| 1603 | h1_ClusterDisp->Fill(aodCluster->GetDispersion()); |
| 1604 | h2_E_time ->Fill(aodCluster->E(),aodCluster->GetTOF()); |
| 1605 | |
| 1606 | // ################################################# |
| 1607 | // track matching eop histograms are handled here... |
| 1608 | // ################################################# |
| 1609 | |
| 1610 | //uncomment this to do the track matching (3 of 3 lines, aod part)!! |
| 1611 | //if(isGoodAodCluster(aodCluster) && !MatchesToTrack){ |
| 1612 | if(isGoodAodCluster(aodCluster)){ |
| 1613 | recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator); |
| 1614 | E1 = aodCluster->E()*recalScale;// TOTAL HACK - JJ |
| 1615 | fAodEv->GetVertex(0)->GetXYZ(vertex); |
| 1616 | aodCluster->GetMomentum(Photon1,vertex); |
| 1617 | Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ |
| 1618 | Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ |
| 1619 | Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ |
| 1620 | Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) ); |
| 1621 | h1_E->Fill(E1); |
| 1622 | h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta()); |
| 1623 | h2_PhiEtaMaxCellCut->Fill(cellphi,celleta); |
| 1624 | } |
| 1625 | clusterPosition.Delete(); |
| 1626 | }//if(aodCluster->IsEMCAL()) |
| 1627 | }//if(fAodEv) |
| 1628 | |
| 1629 | }//loop over nclusters. |
| 1630 | |
| 1631 | //Make same event pions... |
| 1632 | for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){ |
| 1633 | for(unsigned int j=i+1; j<Photons[0][izvtx][imult].size(); j++){ |
| 1634 | Parent = Photons[0][izvtx][imult][i] + Photons[0][izvtx][imult][j]; |
| 1635 | Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]); |
| 1636 | Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]); |
| 1637 | Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[0][izvtx][imult][j].Pt())/ |
| 1638 | (Photons[0][izvtx][imult][i].Pt()+Photons[0][izvtx][imult][j].Pt()); |
| 1639 | Int_t asymCut = 0; |
| 1640 | if (pairasym<0.1) asymCut = 1; |
| 1641 | else if(pairasym<0.7) asymCut = 2; |
| 1642 | else asymCut = 3; |
| 1643 | |
| 1644 | h1_M ->Fill(Parent.M()); |
| 1645 | h3_MptAsymm ->Fill(Parent.M(),Parent.Pt(),asymCut); |
| 1646 | h2_dphi_deta->Fill(deltaphi,deltaeta); |
| 1647 | } |
| 1648 | } |
| 1649 | |
| 1650 | //Make mixed event... |
| 1651 | for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){ |
| 1652 | for(unsigned int ipool=1; ipool<poolDepth; ipool++){ |
| 1653 | for(unsigned int j=0; j<Photons[ipool][izvtx][imult].size(); j++){ |
| 1654 | iskip = randy.Integer(Photons[0][izvtx][imult].size()); |
| 1655 | if(j==iskip) continue; |
| 1656 | Parent = Photons[0][izvtx][imult][i]+Photons[ipool][izvtx][imult][j]; |
| 1657 | Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]); |
| 1658 | Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]); |
| 1659 | Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[ipool][izvtx][imult][j].Pt())/ |
| 1660 | (Photons[0][izvtx][imult][i].Pt()+Photons[ipool][izvtx][imult][j].Pt()); |
| 1661 | Int_t asymCut = 0; |
| 1662 | if (pairasym<0.1) asymCut = 1; |
| 1663 | else if(pairasym<0.7) asymCut = 2; |
| 1664 | else asymCut = 3; |
| 1665 | |
| 1666 | h1_M_mix ->Fill(Parent.M()); |
| 1667 | h3_MptAsymm_mix ->Fill(Parent.M(),Parent.Pt(),asymCut); |
| 1668 | h2_dphi_deta_mix->Fill(deltaphi,deltaeta); |
| 1669 | } |
| 1670 | } |
| 1671 | } |
| 1672 | |
| 1673 | for(int ipool=poolDepth-1; ipool>0; ipool--){ |
| 1674 | Photons[ipool][izvtx][imult].clear(); |
| 1675 | for(unsigned int i=0; i<Photons[ipool-1][izvtx][imult].size(); i++) |
| 1676 | Photons[ipool][izvtx][imult].push_back(Photons[ipool-1][izvtx][imult][i]); |
| 1677 | } |
| 1678 | Photons[0][izvtx][imult].clear(); |
| 1679 | |
| 1680 | |
| 1681 | |
| 1682 | // NEW HISTO should be filled before this point, as PostData puts the |
| 1683 | // information for this iteration of the UserExec in the container |
| 1684 | PostData(1, fOutput); |
| 1685 | } |
| 1686 | |
| 1687 | //________________________________________________________________________ |
| 1688 | void AliAnalysisTaskEMCALMesonGGSDM::Terminate(Option_t *) //specify what you want to have done |
| 1689 | { |
| 1690 | // Called once at the end of the query. |
| 1691 | |
| 1692 | } |
| 1693 | |
| 1694 | //________________________________________________________________________ |
| 1695 | Int_t AliAnalysisTaskEMCALMesonGGSDM::GetZvtxBin(Double_t vertZ) |
| 1696 | { |
| 1697 | |
| 1698 | int izvtx = -1; |
| 1699 | |
| 1700 | if (vertZ<-35) |
| 1701 | izvtx=0; |
| 1702 | else if(vertZ<-30) |
| 1703 | izvtx=1; |
| 1704 | else if(vertZ<-25) |
| 1705 | izvtx=2; |
| 1706 | else if(vertZ<-20) |
| 1707 | izvtx=3; |
| 1708 | else if(vertZ<-15) |
| 1709 | izvtx=4; |
| 1710 | else if(vertZ<-10) |
| 1711 | izvtx=5; |
| 1712 | else if(vertZ< -5) |
| 1713 | izvtx=6; |
| 1714 | else if(vertZ< 0) |
| 1715 | izvtx=7; |
| 1716 | else if(vertZ< 5) |
| 1717 | izvtx=8; |
| 1718 | else if(vertZ< 10) |
| 1719 | izvtx=9; |
| 1720 | else if(vertZ< 15) |
| 1721 | izvtx=10; |
| 1722 | else if(vertZ< 20) |
| 1723 | izvtx=11; |
| 1724 | else if(vertZ< 25) |
| 1725 | izvtx=12; |
| 1726 | else if(vertZ< 30) |
| 1727 | izvtx=13; |
| 1728 | else if(vertZ< 35) |
| 1729 | izvtx=14; |
| 1730 | else |
| 1731 | izvtx=15; |
| 1732 | |
| 1733 | return izvtx; |
| 1734 | } |
| 1735 | |
| 1736 | //________________________________________________________________________ |
| 1737 | Int_t AliAnalysisTaskEMCALMesonGGSDM::GetMultBin(Int_t mult){ |
| 1738 | |
| 1739 | int imult = -1; |
| 1740 | |
| 1741 | if (mult<2) |
| 1742 | imult=0; |
| 1743 | else if(mult<25) |
| 1744 | imult=mult-2; |
| 1745 | else |
| 1746 | imult=24; |
| 1747 | |
| 1748 | return imult; |
| 1749 | } |
| 1750 | |
| 1751 | //________________________________________________________________________ |
| 1752 | Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodEsdCluster(AliESDCaloCluster* esdclust){ |
| 1753 | |
| 1754 | int pass = 1; |
| 1755 | int nMinCells = 2; |
| 1756 | double MinE = 0.4; |
| 1757 | //double MinErat = 0; |
| 1758 | //double MinEcc = 0; |
| 1759 | |
| 1760 | if (!esdclust) |
| 1761 | pass = 0; |
| 1762 | if (!esdclust->IsEMCAL()) |
| 1763 | pass = 0; |
| 1764 | if (esdclust->E()<MinE) |
| 1765 | pass = 0; |
| 1766 | if (esdclust->GetNCells()<nMinCells) |
| 1767 | pass = 0; |
| 1768 | //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat) |
| 1769 | //pass = 0; |
| 1770 | //if (esdclust->Chi2()<MinEcc) // eccentricity cut |
| 1771 | //pass = 0;//this is always -1. |
| 1772 | |
| 1773 | //if(esdclust->GetM02()<0.1) |
| 1774 | // pass = 0; |
| 1775 | //if(esdclust->GetM02()>0.5) |
| 1776 | // pass = 0; |
| 1777 | |
| 1778 | Float_t pos[3] = {0,0,0}; |
| 1779 | esdclust->GetPosition(pos); |
| 1780 | TVector3 clusterPosition(pos); |
| 1781 | if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax || |
| 1782 | clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax ) |
| 1783 | pass = 0; |
| 1784 | clusterPosition.Delete(); |
| 1785 | |
| 1786 | //DOING THIS BY HAND NOW... |
| 1787 | //if(!esdclust->GetNTracksMatched()==0) |
| 1788 | //pass = 0; |
| 1789 | |
| 1790 | return pass; |
| 1791 | } |
| 1792 | |
| 1793 | //________________________________________________________________________ |
| 1794 | Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodAodCluster(AliAODCaloCluster* aodclust){ |
| 1795 | |
| 1796 | int pass = 1; |
| 1797 | int nMinCells = 2; |
| 1798 | double MinE = 0.4; |
| 1799 | //double MinErat = 0; |
| 1800 | //double MinEcc = 0; |
| 1801 | |
| 1802 | if (!aodclust) |
| 1803 | pass = 0; |
| 1804 | if (!aodclust->IsEMCAL()) |
| 1805 | pass = 0; |
| 1806 | if (aodclust->E()<MinE) |
| 1807 | pass = 0; |
| 1808 | if (aodclust->GetNCells()<nMinCells) |
| 1809 | pass = 0; |
| 1810 | //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat) |
| 1811 | //pass = 0; |
| 1812 | //if (aodclust->Chi2()<MinEcc) // eccentricity cut |
| 1813 | //pass = 0;//this is always -1. |
| 1814 | |
| 1815 | //if(aodclust->GetM02()<0.1) |
| 1816 | //pass = 0; |
| 1817 | //if(aodclust->GetM02()>0.5) |
| 1818 | //pass = 0; |
| 1819 | |
| 1820 | Float_t pos[3] = {0,0,0}; |
| 1821 | aodclust->GetPosition(pos); |
| 1822 | TVector3 clusterPosition(pos); |
| 1823 | if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax || |
| 1824 | clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax ) |
| 1825 | pass = 0; |
| 1826 | clusterPosition.Delete(); |
| 1827 | |
| 1828 | //DOING THIS BY HAND NOW... |
| 1829 | //if(!aodclust->GetNTracksMatched()==0) |
| 1830 | //pass = 0; |
| 1831 | |
| 1832 | return pass; |
| 1833 | } |
| 1834 | |
| 1835 | //________________________________________________________________________ |
| 1836 | Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){ |
| 1837 | |
| 1838 | double dphi = p1.Phi() - p2.Phi(); |
| 1839 | |
| 1840 | if(dphi<0.5*TMath::Pi()) |
| 1841 | dphi = dphi + 2.0*TMath::Pi(); |
| 1842 | |
| 1843 | if(dphi>1.5*TMath::Pi()) |
| 1844 | dphi = dphi - 2.0*TMath::Pi(); |
| 1845 | |
| 1846 | return dphi; |
| 1847 | } |
| 1848 | |
| 1849 | //________________________________________________________________________ |
| 1850 | Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaEta(TLorentzVector p1, TLorentzVector p2){ |
| 1851 | |
| 1852 | double deta = p1.PseudoRapidity() - p2.PseudoRapidity(); |
| 1853 | |
| 1854 | return deta; |
| 1855 | } |
| 1856 | |
| 1857 | |
| 1858 | //________________________________________________________________________ |
| 1859 | Double_t AliAnalysisTaskEMCALMesonGGSDM::PrivateEnergyRecal(Double_t energy, Int_t iCalib){ |
| 1860 | |
| 1861 | double recalibfactor = 0.0; |
| 1862 | |
| 1863 | if(iCalib==0){// no recalibration! |
| 1864 | recalibfactor = 1.0; |
| 1865 | } |
| 1866 | else if(iCalib==1){// just a scale factor: |
| 1867 | recalibfactor = 0.984; |
| 1868 | } |
| 1869 | else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC. |
| 1870 | Double_t p[3] = {0.96968, -2.68720, -0.831607}; |
| 1871 | recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0); |
| 1872 | } |
| 1873 | else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??) |
| 1874 | Double_t p[7] = {1.00000e+00, 3.04925e-02, 4.69043e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.00046e+00}; |
| 1875 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); |
| 1876 | } |
| 1877 | else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??) |
| 1878 | Double_t p[7] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01}; |
| 1879 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); |
| 1880 | } |
| 1881 | else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??) |
| 1882 | Double_t p[10] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01, 0.96968, -2.68720, -0.831607}; |
| 1883 | recalibfactor = ( (p[6]/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) ) / ( p[7] + exp(p[8] + p[9]*energy/2.0) ); |
| 1884 | } |
| 1885 | else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT! |
| 1886 | Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968}; |
| 1887 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); |
| 1888 | } |
| 1889 | else if(iCalib==7){// kPi0MCv3 - in AliROOT! |
| 1890 | Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0}; |
| 1891 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); |
| 1892 | } |
| 1893 | else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??) |
| 1894 | Double_t p[10] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286, 0.964, -3.132, -0.435}; |
| 1895 | //Double_t p[10] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286, 0.96968, -2.68720, -0.831607};//same SDM piece as iCalib==2 |
| 1896 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) * (p[7] + exp(p[8]+p[9]*energy*2.0)); |
| 1897 | } |
| 1898 | else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??) |
| 1899 | Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286}; |
| 1900 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); |
| 1901 | } |
| 1902 | else if(iCalib==10){// Jason played with test beam data |
| 1903 | Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054}; |
| 1904 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); |
| 1905 | } |
| 1906 | else if(iCalib==11){// Jason played with test beam MC |
| 1907 | Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437}; |
| 1908 | recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))); |
| 1909 | } |
| 1910 | |
| 1911 | return recalibfactor; |
| 1912 | } |
| 1913 | |
| 1914 | |
| 1915 | //________________________________________________________________________ |
| 1916 | Double_t AliAnalysisTaskEMCALMesonGGSDM::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const |
| 1917 | { |
| 1918 | // Get maximum energy of attached cell. |
| 1919 | |
| 1920 | id = -1; |
| 1921 | AliVCaloCells *fVCells=NULL; |
| 1922 | if(fEsdEv) fVCells = fEsdEv->GetEMCALCells(); |
| 1923 | else if(fAodEv) fVCells = fAodEv->GetEMCALCells(); |
| 1924 | if(!fVCells) |
| 1925 | return 0; |
| 1926 | |
| 1927 | Double_t maxe = 0; |
| 1928 | Int_t ncells = cluster->GetNCells(); |
| 1929 | for (Int_t i=0; i<ncells; i++) { |
| 1930 | Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i))); |
| 1931 | if (e>maxe) { |
| 1932 | maxe = e; |
| 1933 | id = cluster->GetCellAbsId(i); |
| 1934 | } |
| 1935 | } |
| 1936 | return maxe; |
| 1937 | } |
| 1938 | |
| 1939 | |
| 1940 | //________________________________________________________________________ |
| 1941 | Int_t AliAnalysisTaskEMCALMesonGGSDM::IsPhysPrimJ(AliMCEvent *mcEvent, Int_t iTrack){ |
| 1942 | |
| 1943 | AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack)); |
| 1944 | |
| 1945 | Int_t nPTracks= mcEvent->GetNumberOfPrimaries(); |
| 1946 | |
| 1947 | Int_t isPhysPrimary = 1; |
| 1948 | Int_t ismHF = 0; |
| 1949 | Int_t ismLongLivedOrK = 0; |
| 1950 | |
| 1951 | if(mcP->GetMother()<0)//if it has no mother... |
| 1952 | return isPhysPrimary; |
| 1953 | |
| 1954 | Int_t imTrack = mcP->GetMother(); |
| 1955 | AliMCParticle *mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack)); |
| 1956 | |
| 1957 | if( TMath::Abs(mcPm->PdgCode())<10 )//if mother is a single quark... |
| 1958 | return isPhysPrimary; |
| 1959 | |
| 1960 | |
| 1961 | //############################################ |
| 1962 | //get the PDG digits.... |
| 1963 | int num = mcPm->PdgCode(); |
| 1964 | int RevDigits[10] = {0}; |
| 1965 | int nDigits = 0; |
| 1966 | while (num >= 1){ |
| 1967 | RevDigits[nDigits++] = num%10; |
| 1968 | num = num / 10; |
| 1969 | } |
| 1970 | //############################################## |
| 1971 | |
| 1972 | |
| 1973 | if(RevDigits[3]>3)//Baryons |
| 1974 | ismHF = 1; |
| 1975 | else if(RevDigits[2]>3)//Mesons |
| 1976 | ismHF = 1; |
| 1977 | |
| 1978 | ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode()); |
| 1979 | |
| 1980 | if(!ismHF && ismLongLivedOrK) |
| 1981 | isPhysPrimary = 0; |
| 1982 | else{ // check grandmother, greatgrandmothers, etc... |
| 1983 | while(imTrack >= nPTracks){ |
| 1984 | |
| 1985 | if(mcPm->GetMother()<0)//if it has no mother... |
| 1986 | break; |
| 1987 | |
| 1988 | if( TMath::Abs(mcPm->PdgCode()<10) )//if mother is a single quark... |
| 1989 | return isPhysPrimary; |
| 1990 | |
| 1991 | imTrack = mcPm->GetMother(); |
| 1992 | mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack)); |
| 1993 | |
| 1994 | //############################################ |
| 1995 | //get the PDG digits.... |
| 1996 | num = mcPm->PdgCode(); |
| 1997 | for(int i=0; i<10; i++) RevDigits[i] = 0; |
| 1998 | nDigits = 0; |
| 1999 | while (num >= 1){ |
| 2000 | RevDigits[nDigits++] = num%10; |
| 2001 | num = num / 10; |
| 2002 | } |
| 2003 | //############################################## |
| 2004 | if(RevDigits[3]>3)//Baryons |
| 2005 | ismHF = 1; |
| 2006 | else if(RevDigits[2]>3)//Mesons |
| 2007 | ismHF = 1; |
| 2008 | |
| 2009 | ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode()); |
| 2010 | |
| 2011 | if(!ismHF && ismLongLivedOrK) |
| 2012 | isPhysPrimary = 0; |
| 2013 | |
| 2014 | }//while( >=nPTracks) |
| 2015 | } |
| 2016 | |
| 2017 | return isPhysPrimary; |
| 2018 | } |
| 2019 | |
| 2020 | |
| 2021 | //________________________________________________________________________ |
| 2022 | Int_t AliAnalysisTaskEMCALMesonGGSDM::IsLongLivedOrK(Int_t MyPDGcode){ |
| 2023 | |
| 2024 | Int_t MyFlag = 0; |
| 2025 | |
| 2026 | if( |
| 2027 | (TMath::Abs(MyPDGcode) == 22 ) || // Photon |
| 2028 | (TMath::Abs(MyPDGcode) == 11 ) || // Electron |
| 2029 | (TMath::Abs(MyPDGcode) == 13 ) || // Muon(-) |
| 2030 | (TMath::Abs(MyPDGcode) == 211 ) || // Pion |
| 2031 | (TMath::Abs(MyPDGcode) == 321 ) || // Kaon |
| 2032 | (TMath::Abs(MyPDGcode) == 310 ) || // K0s |
| 2033 | (TMath::Abs(MyPDGcode) == 130 ) || // K0l |
| 2034 | (TMath::Abs(MyPDGcode) == 2212) || // Proton |
| 2035 | (TMath::Abs(MyPDGcode) == 2112) || // Neutron |
| 2036 | (TMath::Abs(MyPDGcode) == 3122) || // Lambda_0 |
| 2037 | (TMath::Abs(MyPDGcode) == 3112) || // Sigma Minus |
| 2038 | (TMath::Abs(MyPDGcode) == 3222) || // Sigma Plus |
| 2039 | (TMath::Abs(MyPDGcode) == 3312) || // Xsi Minus |
| 2040 | (TMath::Abs(MyPDGcode) == 3322) || // Xsi |
| 2041 | (TMath::Abs(MyPDGcode) == 3334) || // Omega |
| 2042 | (TMath::Abs(MyPDGcode) == 12 ) || // Electron Neutrino |
| 2043 | (TMath::Abs(MyPDGcode) == 14 ) || // Muon Neutrino |
| 2044 | (TMath::Abs(MyPDGcode) == 16 ) ) // Tau Neutrino |
| 2045 | MyFlag = 1; |
| 2046 | |
| 2047 | return MyFlag; |
| 2048 | } |
| 2049 | |
| 2050 | |
| 2051 | //________________________________________________________________________ |
| 2052 | Int_t AliAnalysisTaskEMCALMesonGGSDM::IsMyMCHeaderType(Int_t iTrack, char *MyType, AliMCEvent *mcEvent) const |
| 2053 | { |
| 2054 | |
| 2055 | Int_t isMyType = 0; |
| 2056 | |
| 2057 | AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(mcEvent->GenEventHeader()); |
| 2058 | if(!cocktail) |
| 2059 | return 0; |
| 2060 | |
| 2061 | TList *genHeaders = cocktail->GetHeaders(); |
| 2062 | |
| 2063 | Int_t nGenerators = genHeaders->GetEntries(); |
| 2064 | Int_t indexMyType = -1; |
| 2065 | Int_t startParticle=0; |
| 2066 | |
| 2067 | for(Int_t igen = 0; igen < nGenerators; igen++){ |
| 2068 | AliGenEventHeader* eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ; |
| 2069 | TString name = eventHeader2->GetName(); |
| 2070 | startParticle += eventHeader2->NProduced(); |
| 2071 | //cout << name << endl; |
| 2072 | if (name.Contains(MyType,TString::kIgnoreCase)){ |
| 2073 | indexMyType = igen; |
| 2074 | startParticle -= eventHeader2->NProduced(); |
| 2075 | break; |
| 2076 | } |
| 2077 | } |
| 2078 | |
| 2079 | AliGenEventHeader *addedPi0Header = (AliGenEventHeader*)genHeaders->At(indexMyType); |
| 2080 | Int_t ipi0min = startParticle; |
| 2081 | Int_t ipi0max = ipi0min+addedPi0Header->NProduced()-1; |
| 2082 | if(iTrack >= ipi0min && iTrack <= ipi0max) |
| 2083 | isMyType = 1; |
| 2084 | |
| 2085 | return isMyType; |
| 2086 | } |
| 2087 | |
| 2088 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame