| 0bd4b00f |
1 | #include "AliForwardMCCorrectionsTask.h" |
| 563a673f |
2 | #include "AliTriggerAnalysis.h" |
| 3 | #include "AliPhysicsSelection.h" |
| 4 | #include "AliLog.h" |
| 0bd4b00f |
5 | // #include "AliFMDAnaParameters.h" |
| 6 | #include "AliHeader.h" |
| 7 | #include "AliGenEventHeader.h" |
| 563a673f |
8 | #include "AliESDEvent.h" |
| 9 | #include "AliAODHandler.h" |
| 10 | #include "AliMultiplicity.h" |
| 11 | #include "AliInputEventHandler.h" |
| 0bd4b00f |
12 | #include "AliStack.h" |
| 13 | #include "AliMCEvent.h" |
| 14 | #include "AliFMDStripIndex.h" |
| 563a673f |
15 | #include <TH1.h> |
| 0bd4b00f |
16 | #include <TH3D.h> |
| 563a673f |
17 | #include <TDirectory.h> |
| 18 | #include <TTree.h> |
| 19 | |
| 20 | //==================================================================== |
| 21 | namespace { |
| 22 | const char* GetEventName(Bool_t tr, Bool_t vtx) |
| 23 | { |
| 0bd4b00f |
24 | return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : "")); |
| 563a673f |
25 | } |
| 26 | const char* GetHitsName(UShort_t d, Char_t r) |
| 27 | { |
| 28 | return Form("hitsFMD%d%c", d, r); |
| 29 | } |
| 30 | const char* GetStripsName(UShort_t d, Char_t r) |
| 31 | { |
| 32 | return Form("stripsFMD%d%c", d, r); |
| 33 | } |
| 34 | const char* GetPrimaryName(Char_t r, Bool_t trVtx) |
| 35 | { |
| 36 | return Form("primaries%s%s", |
| 37 | ((r == 'I' || r == 'i') ? "Inner" : "Outer"), |
| 38 | (trVtx ? "TrVtx" : "All")); |
| 39 | } |
| 40 | } |
| 41 | |
| 42 | //==================================================================== |
| 0bd4b00f |
43 | AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask() |
| 563a673f |
44 | : AliAnalysisTaskSE(), |
| 0bd4b00f |
45 | fHEvents(0), |
| 563a673f |
46 | fHEventsTr(0), |
| 47 | fHEventsTrVtx(0), |
| 0bd4b00f |
48 | fHEventsVtx(0), |
| 563a673f |
49 | fHTriggers(0), |
| 0bd4b00f |
50 | fPrimaryInnerAll(0), |
| 51 | fPrimaryOuterAll(0), |
| 52 | fPrimaryInnerTrVtx(0), |
| 53 | fPrimaryOuterTrVtx(0), |
| 54 | fHitsFMD1i(0), |
| 55 | fHitsFMD2i(0), |
| 56 | fHitsFMD2o(0), |
| 57 | fHitsFMD3i(0), |
| 58 | fHitsFMD3o(0), |
| 59 | fStripsFMD1i(0), |
| 60 | fStripsFMD2i(0), |
| 61 | fStripsFMD2o(0), |
| 62 | fStripsFMD3i(0), |
| 63 | fStripsFMD3o(0), |
| 563a673f |
64 | fVtxAxis(), |
| 0bd4b00f |
65 | fEtaAxis(), |
| 66 | fList() |
| 563a673f |
67 | { |
| 68 | } |
| 69 | |
| 70 | //____________________________________________________________________ |
| 0bd4b00f |
71 | AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const char* name) |
| 563a673f |
72 | : AliAnalysisTaskSE(name), |
| 0bd4b00f |
73 | fHEvents(0), |
| 563a673f |
74 | fHEventsTr(0), |
| 0bd4b00f |
75 | fHEventsTrVtx(0), |
| 76 | fHEventsVtx(0), |
| 563a673f |
77 | fHTriggers(0), |
| 0bd4b00f |
78 | fPrimaryInnerAll(0), |
| 79 | fPrimaryOuterAll(0), |
| 80 | fPrimaryInnerTrVtx(0), |
| 81 | fPrimaryOuterTrVtx(0), |
| 82 | fHitsFMD1i(0), |
| 83 | fHitsFMD2i(0), |
| 84 | fHitsFMD2o(0), |
| 85 | fHitsFMD3i(0), |
| 86 | fHitsFMD3o(0), |
| 87 | fStripsFMD1i(0), |
| 88 | fStripsFMD2i(0), |
| 89 | fStripsFMD2o(0), |
| 90 | fStripsFMD3i(0), |
| 91 | fStripsFMD3o(0), |
| 563a673f |
92 | fVtxAxis(10,-10,10), |
| 0bd4b00f |
93 | fEtaAxis(200,-4,6), |
| 94 | fList() |
| 563a673f |
95 | { |
| 96 | DefineOutput(1, TList::Class()); |
| 97 | } |
| 98 | |
| 99 | //____________________________________________________________________ |
| 0bd4b00f |
100 | AliForwardMCCorrectionsTask::AliForwardMCCorrectionsTask(const AliForwardMCCorrectionsTask& o) |
| 563a673f |
101 | : AliAnalysisTaskSE(o), |
| 0bd4b00f |
102 | fHEvents(o.fHEvents), |
| 563a673f |
103 | fHEventsTr(o.fHEventsTr), |
| 0bd4b00f |
104 | fHEventsTrVtx(o.fHEventsTrVtx), |
| 105 | fHEventsVtx(o.fHEventsVtx), |
| 563a673f |
106 | fHTriggers(o.fHTriggers), |
| 0bd4b00f |
107 | fPrimaryInnerAll(o.fPrimaryInnerAll), |
| 108 | fPrimaryOuterAll(o.fPrimaryOuterAll), |
| 109 | fPrimaryInnerTrVtx(o.fPrimaryInnerTrVtx), |
| 110 | fPrimaryOuterTrVtx(o.fPrimaryOuterTrVtx), |
| 111 | fHitsFMD1i(o.fHitsFMD1i), |
| 112 | fHitsFMD2i(o.fHitsFMD2i), |
| 113 | fHitsFMD2o(o.fHitsFMD2o), |
| 114 | fHitsFMD3i(o.fHitsFMD3i), |
| 115 | fHitsFMD3o(o.fHitsFMD3o), |
| 116 | fStripsFMD1i(o.fStripsFMD1i), |
| 117 | fStripsFMD2i(o.fStripsFMD2i), |
| 118 | fStripsFMD2o(o.fStripsFMD2o), |
| 119 | fStripsFMD3i(o.fStripsFMD3i), |
| 120 | fStripsFMD3o(o.fStripsFMD3o), |
| 121 | fVtxAxis(o.fVtxAxis.GetNbins(), o.fVtxAxis.GetXmin(), o.fVtxAxis.GetXmax()), |
| 122 | fEtaAxis(o.fEtaAxis.GetNbins(), o.fEtaAxis.GetXmin(), o.fEtaAxis.GetXmax()), |
| 123 | fList(o.fList) |
| 563a673f |
124 | { |
| 125 | } |
| 126 | |
| 127 | //____________________________________________________________________ |
| 0bd4b00f |
128 | AliForwardMCCorrectionsTask& |
| 129 | AliForwardMCCorrectionsTask::operator=(const AliForwardMCCorrectionsTask& o) |
| 563a673f |
130 | { |
| 131 | fHEventsTr = o.fHEventsTr; |
| 132 | fHEventsTrVtx = o.fHEventsTrVtx; |
| 133 | fHTriggers = o.fHTriggers; |
| 0bd4b00f |
134 | SetVertexAxis(o.fVtxAxis); |
| 135 | SetEtaAxis(o.fEtaAxis); |
| 563a673f |
136 | |
| 137 | return *this; |
| 138 | } |
| 139 | |
| 140 | //____________________________________________________________________ |
| 141 | void |
| 0bd4b00f |
142 | AliForwardMCCorrectionsTask::Init() |
| 563a673f |
143 | { |
| 144 | } |
| 145 | |
| 146 | //____________________________________________________________________ |
| 147 | void |
| 0bd4b00f |
148 | AliForwardMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min, |
| 149 | Double_t max) |
| 563a673f |
150 | { |
| 151 | if (max < min) max = -min; |
| 152 | if (min < max) { |
| 153 | Double_t tmp = min; |
| 154 | min = max; |
| 155 | max = tmp; |
| 156 | } |
| 157 | if (min < -10) |
| 158 | AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min)); |
| 159 | if (max > +10) |
| 160 | AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max)); |
| 161 | fVtxAxis.Set(nBin, min, max); |
| 162 | } |
| 163 | //____________________________________________________________________ |
| 164 | void |
| 0bd4b00f |
165 | AliForwardMCCorrectionsTask::SetVertexAxis(const TAxis& axis) |
| 563a673f |
166 | { |
| 167 | SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax()); |
| 168 | } |
| 169 | |
| 170 | //____________________________________________________________________ |
| 171 | void |
| 0bd4b00f |
172 | AliForwardMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max) |
| 563a673f |
173 | { |
| 174 | if (max < min) max = -min; |
| 175 | if (min < max) { |
| 176 | Double_t tmp = min; |
| 177 | min = max; |
| 178 | max = tmp; |
| 179 | } |
| 180 | if (min < -4) |
| 181 | AliWarning(Form("Minimum eta %f < -4, make sure you want this",min)); |
| 182 | if (max > +6) |
| 183 | AliWarning(Form("Minimum eta %f > +6, make sure you want this",max)); |
| 184 | fVtxAxis.Set(nBin, min, max); |
| 185 | } |
| 186 | //____________________________________________________________________ |
| 187 | void |
| 0bd4b00f |
188 | AliForwardMCCorrectionsTask::SetEtaAxis(const TAxis& axis) |
| 563a673f |
189 | { |
| 190 | SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax()); |
| 191 | } |
| 192 | |
| 563a673f |
193 | //____________________________________________________________________ |
| 194 | TH3D* |
| 0bd4b00f |
195 | AliForwardMCCorrectionsTask::Make3D(const char* name, const char* title, |
| 563a673f |
196 | Int_t nPhi) const |
| 197 | { |
| 198 | TH3D* ret = new TH3D(name, title, |
| 199 | fVtxAxis.GetNbins(), |
| 200 | fVtxAxis.GetXmin(), |
| 201 | fVtxAxis.GetXmax(), |
| 202 | fEtaAxis.GetNbins(), |
| 203 | fEtaAxis.GetXmin(), |
| 204 | fEtaAxis.GetXmax(), |
| 205 | nPhi, 0, 2*TMath::Pi()); |
| 206 | ret->SetXTitle("v_{z} [cm]"); |
| 207 | ret->SetYTitle("#eta"); |
| 208 | ret->SetZTitle("#varphi [radians]"); |
| 209 | ret->SetDirectory(0); |
| 210 | ret->SetStats(0); |
| 211 | ret->Sumw2(); |
| 212 | |
| 213 | return ret; |
| 214 | } |
| 215 | //____________________________________________________________________ |
| 216 | TH1D* |
| 0bd4b00f |
217 | AliForwardMCCorrectionsTask::Make1D(const char* name, const char* title) const |
| 563a673f |
218 | { |
| 219 | TH1D* ret = new TH1D(name, title, |
| 220 | fEtaAxis.GetNbins(), |
| 221 | fEtaAxis.GetXmin(), |
| 222 | fEtaAxis.GetXmax()); |
| 223 | ret->SetXTitle("#eta"); |
| 224 | ret->SetFillColor(kRed+1); |
| 225 | ret->SetFillStyle(3001); |
| 226 | ret->SetDirectory(0); |
| 227 | ret->SetStats(0); |
| 228 | ret->Sumw2(); |
| 229 | |
| 230 | return ret; |
| 231 | } |
| 232 | |
| 233 | //____________________________________________________________________ |
| 234 | void |
| 0bd4b00f |
235 | AliForwardMCCorrectionsTask::UserCreateOutputObjects() |
| 563a673f |
236 | { |
| 237 | fList = new TList; |
| 238 | fList->SetName(GetName()); |
| 239 | |
| 0bd4b00f |
240 | fHEvents = new TH1I(GetEventName(false,false), |
| 563a673f |
241 | "Number of all events", |
| 0bd4b00f |
242 | fVtxAxis.GetNbins(), |
| 243 | fVtxAxis.GetXmin(), |
| 244 | fVtxAxis.GetXmax()); |
| 563a673f |
245 | fHEvents->SetXTitle("v_{z} [cm]"); |
| 246 | fHEvents->SetYTitle("# of events"); |
| 247 | fHEvents->SetFillColor(kBlue+1); |
| 248 | fHEvents->SetFillStyle(3001); |
| 249 | fHEvents->SetDirectory(0); |
| 250 | fList->Add(fHEvents); |
| 251 | |
| 0bd4b00f |
252 | fHEventsTr = new TH1I(GetEventName(true, false), |
| 253 | "Number of triggered events", |
| 563a673f |
254 | fVtxAxis.GetNbins(), |
| 255 | fVtxAxis.GetXmin(), |
| 256 | fVtxAxis.GetXmax()); |
| 257 | fHEventsTr->SetXTitle("v_{z} [cm]"); |
| 258 | fHEventsTr->SetYTitle("# of events"); |
| 259 | fHEventsTr->SetFillColor(kRed+1); |
| 260 | fHEventsTr->SetFillStyle(3001); |
| 261 | fHEventsTr->SetDirectory(0); |
| 262 | fList->Add(fHEventsTr); |
| 263 | |
| 0bd4b00f |
264 | fHEventsTrVtx = new TH1I(GetEventName(true,true), |
| 563a673f |
265 | "Number of events w/trigger and vertex", |
| 266 | fVtxAxis.GetNbins(), |
| 267 | fVtxAxis.GetXmin(), |
| 268 | fVtxAxis.GetXmax()); |
| 269 | fHEventsTrVtx->SetXTitle("v_{z} [cm]"); |
| 270 | fHEventsTrVtx->SetYTitle("# of events"); |
| 271 | fHEventsTrVtx->SetFillColor(kBlue+1); |
| 272 | fHEventsTrVtx->SetFillStyle(3001); |
| 273 | fHEventsTrVtx->SetDirectory(0); |
| 274 | fList->Add(fHEventsTrVtx); |
| 0bd4b00f |
275 | |
| 276 | fHEventsVtx = new TH1I(GetEventName(false,true), |
| 563a673f |
277 | "Number of events w/vertex", |
| 278 | fVtxAxis.GetNbins(), |
| 279 | fVtxAxis.GetXmin(), |
| 280 | fVtxAxis.GetXmax()); |
| 281 | fHEventsVtx->SetXTitle("v_{z} [cm]"); |
| 282 | fHEventsVtx->SetYTitle("# of events"); |
| 283 | fHEventsVtx->SetFillColor(kBlue+1); |
| 284 | fHEventsVtx->SetFillStyle(3001); |
| 285 | fHEventsVtx->SetDirectory(0); |
| 286 | fList->Add(fHEventsVtx); |
| 287 | |
| 288 | |
| 289 | fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10); |
| 290 | fHTriggers->SetFillColor(kRed+1); |
| 291 | fHTriggers->SetFillStyle(3001); |
| 292 | fHTriggers->SetStats(0); |
| 293 | fHTriggers->SetDirectory(0); |
| 294 | fHTriggers->GetXaxis()->SetBinLabel(1,"INEL"); |
| 295 | fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0"); |
| 296 | fHTriggers->GetXaxis()->SetBinLabel(3,"NSD"); |
| 297 | fHTriggers->GetXaxis()->SetBinLabel(4,"Empty"); |
| 298 | fHTriggers->GetXaxis()->SetBinLabel(5,"A"); |
| 299 | fHTriggers->GetXaxis()->SetBinLabel(6,"B"); |
| 300 | fHTriggers->GetXaxis()->SetBinLabel(7,"C"); |
| 301 | fHTriggers->GetXaxis()->SetBinLabel(8,"E"); |
| 302 | fList->Add(fHTriggers); |
| 303 | |
| 304 | fPrimaryInnerAll = Make3D(GetPrimaryName('I',false), "Primary particles, " |
| 305 | "20 #varphi bins, all events", 20); |
| 306 | fPrimaryOuterAll = Make3D(GetPrimaryName('O',false), "Primary particles, " |
| 307 | "40 #varphi bins, all events", 40); |
| 308 | fPrimaryInnerTrVtx = Make3D(GetPrimaryName('I',true), "Primary particles, " |
| 309 | "20 #varphi bins, Tr+Vtx events", 20); |
| 310 | fPrimaryOuterTrVtx = Make3D(GetPrimaryName('O',true), "Primary particles, " |
| 311 | "40 #varphi bins, Tr+Vtx events", 40); |
| 312 | TList* primaries = new TList; |
| 313 | primaries->SetName("primaries"); |
| 314 | primaries->Add(fPrimaryInnerAll); |
| 315 | primaries->Add(fPrimaryOuterAll); |
| 316 | primaries->Add(fPrimaryInnerTrVtx); |
| 317 | primaries->Add(fPrimaryOuterTrVtx); |
| 318 | fList->Add(primaries); |
| 319 | |
| 320 | |
| 321 | fHitsFMD1i = Make3D(GetHitsName(1,'i'), "All hits in FMD1i, Tr+Vtx", 20); |
| 322 | fHitsFMD2i = Make3D(GetHitsName(2,'i'), "All hits in FMD2i, Tr+Vtx", 20); |
| 323 | fHitsFMD2o = Make3D(GetHitsName(2,'o'), "All hits in FMD2o, Tr+Vtx", 40); |
| 324 | fHitsFMD3i = Make3D(GetHitsName(3,'i'), "All hits in FMD3i, Tr+Vtx", 20); |
| 325 | fHitsFMD3o = Make3D(GetHitsName(3,'o'), "All hits in FMD3o, Tr+Vtx", 40); |
| 326 | TList* hits = new TList; |
| 327 | hits->SetName("hits"); |
| 328 | hits->Add(fHitsFMD1i); |
| 329 | hits->Add(fHitsFMD2i); |
| 330 | hits->Add(fHitsFMD2o); |
| 331 | hits->Add(fHitsFMD3i); |
| 332 | hits->Add(fHitsFMD3o); |
| 333 | fList->Add(hits); |
| 334 | |
| 335 | fStripsFMD1i = Make1D(GetStripsName(1,'i'), "# hit strips in FMD1i, Tr+Vtx"); |
| 336 | fStripsFMD2i = Make1D(GetStripsName(2,'i'), "# hit strips in FMD2i, Tr+Vtx"); |
| 337 | fStripsFMD2o = Make1D(GetStripsName(2,'o'), "# hit strips in FMD2o, Tr+Vtx"); |
| 338 | fStripsFMD3i = Make1D(GetStripsName(3,'i'), "# hit strips in FMD3i, Tr+Vtx"); |
| 339 | fStripsFMD3o = Make1D(GetStripsName(3,'o'), "# hit strips in FMD3o, Tr+Vtx"); |
| 340 | TList* strips = new TList; |
| 341 | strips->SetName("strips"); |
| 342 | strips->Add(fStripsFMD1i); |
| 343 | strips->Add(fStripsFMD2i); |
| 344 | strips->Add(fStripsFMD2o); |
| 345 | strips->Add(fStripsFMD3i); |
| 346 | strips->Add(fStripsFMD3o); |
| 347 | fList->Add(strips); |
| 348 | |
| 349 | PostData(1, fList); |
| 350 | } |
| 351 | //____________________________________________________________________ |
| 352 | void |
| 0bd4b00f |
353 | AliForwardMCCorrectionsTask::UserExec(Option_t*) |
| 563a673f |
354 | { |
| 355 | // Get the input data - MC event |
| 356 | AliMCEvent* mcEvent = MCEvent(); |
| 0bd4b00f |
357 | if (mcEvent) { |
| 563a673f |
358 | AliWarning("No MC event found"); |
| 359 | return; |
| 360 | } |
| 361 | |
| 362 | // Get the input data - ESD event |
| 363 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()); |
| 364 | if (!esd) { |
| 365 | AliWarning("No ESD event found for input event"); |
| 366 | return; |
| 367 | } |
| 368 | |
| 369 | // Get the particle stack |
| 370 | AliStack* stack = mcEvent->Stack(); |
| 371 | |
| 372 | // Get the event generate header |
| 373 | AliHeader* mcHeader = mcEvent->Header(); |
| 0bd4b00f |
374 | AliGenEventHeader* genHeader = mcHeader->GenEventHeader(); |
| 563a673f |
375 | |
| 376 | // Get the generator vertex |
| 377 | TArrayF mcVertex; |
| 0bd4b00f |
378 | genHeader->PrimaryVertex(mcVertex); |
| 563a673f |
379 | Double_t mcVtxZ = mcVertex.At(2); |
| 380 | |
| 381 | // Check the MC vertex is in range |
| 0bd4b00f |
382 | Int_t mcVtxBin = fVtxAxis.FindBin(mcVtxZ); |
| 383 | if (mcVtxBin < 1 || mcVtxBin > fVtxAxis.GetNbins()) { |
| 563a673f |
384 | #ifdef VERBOSE |
| 385 | AliWarning(Form("MC v_z=%f is out of rante [%,%f]", |
| 0bd4b00f |
386 | mcVtxBin, fVtxAxis.GetXmin(), fVtxAxis.GetXmax())); |
| 563a673f |
387 | #endif |
| 388 | return; |
| 389 | } |
| 390 | |
| 0bd4b00f |
391 | UInt_t triggers; |
| 392 | Bool_t gotTrigggers; |
| 393 | Bool_t gotInel; |
| 394 | Double_t vZ; |
| 395 | Bool_t gotVertex; |
| 396 | #if 0 |
| 397 | // Use event inspector instead |
| 563a673f |
398 | // Get the triggers |
| 399 | UInt_t triggers = 0; |
| 400 | Bool_t gotTriggers = AliForwardUtil::ReadTriggers(esd,triggers,fHTriggers); |
| 401 | Bool_t gotInel = triggers & AliAODForwardMult::kInel; |
| 402 | |
| 403 | // Get the ESD vertex |
| 404 | Double_t vZ = -1000000; |
| 405 | Bool_t gotVertex = AliForwardUtil::ReadVertex(esd,vZ); |
| 0bd4b00f |
406 | #endif |
| 407 | |
| 563a673f |
408 | |
| 409 | // Fill the event count histograms |
| 410 | if (gotInel) fHEventsTr->Fill(mcVtxZ); |
| 411 | if (gotInel && gotVertex) fHEventsTrVtx->Fill(mcVtxZ); |
| 412 | if (gotVertex) fHEventsVtx->Fill(mcVtxZ); |
| 413 | fHEvents->Fill(mcVtxZ); |
| 414 | |
| 415 | // Cache of the hits |
| 416 | AliFMDFloatMap hitsByStrip; |
| 417 | AliFMDFloatMap lastByStrip; |
| 418 | |
| 419 | // Loop over all tracks |
| 420 | Int_t nTracks = mcEvent->GetNumberOfTracks(); |
| 421 | for (Int_t iTr = 0; iTr < nTracks; iTr++) { |
| 422 | AliMCParticle* particle = |
| 423 | static_cast<AliMCParticle*>(mcEvent->GetTrack(iTr)); |
| 424 | |
| 425 | // Check the returned particle |
| 426 | if (!particle) continue; |
| 427 | |
| 428 | // Check if this charged and a primary |
| 429 | Bool_t isCharged = particle->Charge() != 0; |
| 430 | Bool_t isPrimary = stack->IsPhysicalPrimary(iTr); |
| 431 | |
| 432 | // Fill (eta,phi) of the particle into histograsm for b |
| 433 | Double_t eta = particle->Eta(); |
| 0bd4b00f |
434 | Double_t phi = particle->Phi(); |
| 563a673f |
435 | |
| 436 | if (isCharged && isPrimary) |
| 437 | FillPrimary(gotInel, gotVertex, mcVtxZ, eta, phi); |
| 438 | |
| 439 | // For the rest, ignore non-collisions, and events out of vtx range |
| 0bd4b00f |
440 | if (!gotInel || gotVertex) continue; |
| 563a673f |
441 | |
| 442 | Int_t nTrRef = particle->GetNumberOfTrackReferences(); |
| 443 | for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { |
| 444 | AliTrackReference* trackRef = particle->GetTrackReference(iTrRef); |
| 445 | |
| 446 | // Check existence |
| 0bd4b00f |
447 | if (!trackRef) continue; |
| 563a673f |
448 | |
| 449 | // Check that we hit an FMD element |
| 0bd4b00f |
450 | if (trackRef->DetectorId() != AliTrackReference::kFMD) |
| 563a673f |
451 | continue; |
| 452 | |
| 453 | // Get the detector coordinates |
| 454 | UShort_t d, s, t; |
| 455 | Char_t r; |
| 0bd4b00f |
456 | AliFMDStripIndex::Unpack(trackRef->UserId(), d, r, s, t); |
| 563a673f |
457 | |
| 458 | // Check if mother (?) is charged and that we haven't dealt with it |
| 459 | // already |
| 460 | Int_t lastTrack = Int_t(lastByStrip(d,r,s,t)); |
| 461 | if (!isCharged || iTr == lastTrack) continue; |
| 462 | |
| 463 | // Increase counter of hits per strip |
| 464 | hitsByStrip(d,r,s,t) += 1; |
| 465 | |
| 0bd4b00f |
466 | // Double_t trRefEta = esd->GetFMDData()->Eta(d,r,s,t); |
| 467 | // Double_t trRefPhi = esd->GetFMDData()->Phi(d,r,s,t); |
| 563a673f |
468 | |
| 469 | // Fill strip information into histograms |
| 0bd4b00f |
470 | FillStrip(d, r, mcVtxZ, eta, phi, hitsByStrip(d,r,s,t) == 1); |
| 563a673f |
471 | |
| 472 | // Set the last processed track number - marking it as done for |
| 473 | // this strip |
| 474 | lastByStrip(d,r,s,t) = Float_t(iTr); |
| 475 | |
| 476 | // Flag neighboring strips too |
| 477 | Int_t nStrip = (r == 'I' || r == 'i' ? 512 : 256); |
| 478 | if (t > 0) lastByStrip(d,r,s,t-1) = Float_t(iTr); |
| 479 | if (t < nStrip-1) lastByStrip(d,r,s,t+1) = Float_t(iTr); |
| 480 | } |
| 481 | } |
| 482 | } |
| 483 | |
| 484 | //____________________________________________________________________ |
| 485 | void |
| 0bd4b00f |
486 | AliForwardMCCorrectionsTask::FillPrimary(Bool_t gotInel, Bool_t gotVtx, |
| 563a673f |
487 | Double_t vz, Double_t eta, Double_t phi) |
| 488 | { |
| 489 | if (gotInel && gotVtx) { |
| 490 | // This takes the place of hPrimary_FMD_<r>_vtx<v> and |
| 491 | // Analysed_FMD<r>_vtx<v> |
| 492 | fPrimaryInnerTrVtx->Fill(vz,eta,phi); |
| 493 | fPrimaryOuterTrVtx->Fill(vz,eta,phi); |
| 494 | } |
| 495 | // This takes the place of Inel_FMD<r>_vtx<v> and |
| 496 | // Analysed_FMD<r>_vtx<v> |
| 497 | fPrimaryInnerAll->Fill(vz,eta,phi); |
| 498 | fPrimaryOuterAll->Fill(vz,eta,phi); |
| 499 | } |
| 500 | |
| 501 | //____________________________________________________________________ |
| 502 | void |
| 0bd4b00f |
503 | AliForwardMCCorrectionsTask::FillStrip(UShort_t d, Char_t r, |
| 563a673f |
504 | Double_t vz, Double_t eta, Double_t phi, |
| 505 | Bool_t first) |
| 506 | { |
| 507 | // Number of hit strips per eta bin |
| 0bd4b00f |
508 | TH1D* strips = 0; |
| 563a673f |
509 | // All hits per ring, vertex in eta,phi bins. This takes the place |
| 510 | // of hHits_FMD<d><r>_vtx<v> and DoubleHits_FMD<d><r> (the later |
| 511 | // being just the 2D projection over the X bins) |
| 0bd4b00f |
512 | TH3D* hits = 0; |
| 563a673f |
513 | switch (d) { |
| 514 | case 1: |
| 515 | hits = fHitsFMD1i; |
| 516 | strips = fStripsFMD1i; |
| 517 | break; |
| 518 | case 2: |
| 519 | hits = (r == 'I' || r == 'i' ? fHitsFMD2i : fHitsFMD2o); |
| 520 | strips = (r == 'I' || r == 'i' ? fStripsFMD2i : fStripsFMD2o); |
| 521 | break; |
| 522 | case 3: |
| 523 | hits = (r == 'I' || r == 'i' ? fHitsFMD3i : fHitsFMD3o); |
| 524 | strips = (r == 'I' || r == 'i' ? fStripsFMD3i : fStripsFMD3o); |
| 525 | break; |
| 526 | } |
| 527 | if (!hits || !strips) return; |
| 528 | |
| 529 | if (first) strips->Fill(eta); |
| 530 | |
| 531 | hits->Fill(vz, eta, phi); |
| 532 | } |
| 533 | |
| 534 | //____________________________________________________________________ |
| 535 | TH2D* |
| 0bd4b00f |
536 | AliForwardMCCorrectionsTask::GetVertexProj(Int_t v, TH3D* src) const |
| 563a673f |
537 | { |
| 538 | Int_t nX = src->GetNbinsX(); |
| 539 | if (v > nX || v < 1) return 0; |
| 540 | |
| 541 | src->GetXaxis()->SetRange(v,v+1); |
| 542 | |
| 543 | TH2D* ret = static_cast<TH2D*>(src->Project3D("zy")); |
| 544 | ret->SetName(Form("%s_vtx%02d", src->GetName(), v)); |
| 545 | ret->SetDirectory(0); |
| 546 | |
| 547 | src->GetXaxis()->SetRange(0,nX+1); |
| 548 | |
| 549 | return ret; |
| 550 | } |
| 551 | |
| 552 | |
| 553 | //____________________________________________________________________ |
| 554 | void |
| 0bd4b00f |
555 | AliForwardMCCorrectionsTask::Terminate(Option_t*) |
| 563a673f |
556 | { |
| 557 | TList* list = dynamic_cast<TList*>(GetOutputData(1)); |
| 558 | if (!list) { |
| 559 | AliError("No output list defined"); |
| 560 | return; |
| 561 | } |
| 562 | |
| 563 | TList* primaries = static_cast<TList*>(list->FindObject("primaries")); |
| 564 | TList* hits = static_cast<TList*>(list->FindObject("hits")); |
| 565 | TList* strips = static_cast<TList*>(list->FindObject("strips")); |
| 566 | if (!primaries || !hits || !strips) { |
| 567 | AliError(Form("Could not find all sub-lists in %s (p=%p,h=%p,s=%p)", |
| 568 | list->GetName(), primaries, hits, strips)); |
| 569 | return; |
| 570 | } |
| 571 | |
| 572 | TH1I* eventsAll = |
| 0bd4b00f |
573 | static_cast<TH1I*>(list->FindObject(GetEventName(false,false))); |
| 563a673f |
574 | TH1I* eventsTr = |
| 0bd4b00f |
575 | static_cast<TH1I*>(list->FindObject(GetEventName(true,false))); |
| 563a673f |
576 | TH1I* eventsVtx = |
| 0bd4b00f |
577 | static_cast<TH1I*>(list->FindObject(GetEventName(false,true))); |
| 563a673f |
578 | TH1I* eventsTrVtx = |
| 0bd4b00f |
579 | static_cast<TH1I*>(list->FindObject(GetEventName(true,true))); |
| 563a673f |
580 | if (!eventsAll || !eventsTr || !eventsVtx || !eventsTrVtx) { |
| 581 | AliError(Form("Could not find all event histograms in %s " |
| 582 | "(a=%p,t=%p,v=%p,tv=%p)", list->GetName(), |
| 583 | eventsAll, eventsTr, eventsVtx, eventsTrVtx)); |
| 584 | return; |
| 585 | } |
| 586 | |
| 587 | TH3D* primIAll = |
| 588 | static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('I',false))); |
| 589 | TH3D* primOAll = |
| 590 | static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',false))); |
| 591 | TH3D* primITrVtx = |
| 592 | static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('I',true))); |
| 593 | TH3D* primOTrVtx = |
| 594 | static_cast<TH3D*>(primaries->FindObject(GetPrimaryName('O',true))); |
| 595 | if (!primIAll || !primOAll || !primITrVtx || !primOTrVtx) { |
| 596 | AliError(Form("Could not find all primary particle histograms in %s " |
| 597 | "(ai=%p,ao=%p,tvi=%p,tvo=%p)", list->GetName(), |
| 598 | primIAll, primOAll, primITrVtx, primOTrVtx)); |
| 599 | return; |
| 600 | } |
| 601 | |
| 602 | TH3D* hits1i = static_cast<TH3D*>(hits->FindObject(GetHitsName(1,'i'))); |
| 603 | TH3D* hits2i = static_cast<TH3D*>(hits->FindObject(GetHitsName(2,'i'))); |
| 604 | TH3D* hits2o = static_cast<TH3D*>(hits->FindObject(GetHitsName(2,'o'))); |
| 605 | TH3D* hits3i = static_cast<TH3D*>(hits->FindObject(GetHitsName(3,'i'))); |
| 606 | TH3D* hits3o = static_cast<TH3D*>(hits->FindObject(GetHitsName(3,'o'))); |
| 607 | if (!hits1i || !hits2i || !hits2o || hits3i || hits3o) { |
| 608 | AliError(Form("Could not find all ring hits histograms in %s " |
| 609 | "(1i=%p,2i=%p,2o=%p,3i=%p,3o=%p)", hits->GetName(), |
| 610 | hits1i, hits2i, hits2o, hits3i, hits3o)); |
| 611 | return; |
| 612 | } |
| 613 | |
| 614 | TH1D* strips1i = static_cast<TH1D*>(strips->FindObject(GetStripsName(1,'i'))); |
| 615 | TH1D* strips2i = static_cast<TH1D*>(strips->FindObject(GetStripsName(2,'i'))); |
| 616 | TH1D* strips2o = static_cast<TH1D*>(strips->FindObject(GetStripsName(2,'o'))); |
| 617 | TH1D* strips3i = static_cast<TH1D*>(strips->FindObject(GetStripsName(3,'i'))); |
| 618 | TH1D* strips3o = static_cast<TH1D*>(strips->FindObject(GetStripsName(3,'o'))); |
| 619 | if (!strips1i || !strips2i || !strips2o || strips3i || strips3o) { |
| 620 | AliError(Form("Could not find all ring strips histograms in %s " |
| 621 | "(1i=%p,2i=%p,2o=%p,3i=%p,3o=%p)", strips->GetName(), |
| 622 | strips1i, strips2i, strips2o, strips3i, strips3o)); |
| 623 | return; |
| 624 | } |
| 625 | |
| 626 | // Calculate the over-all event selection efficiency |
| 627 | TH1D* selEff = new TH1D("selEff", "Event selection efficiency", |
| 0bd4b00f |
628 | fVtxAxis.GetNbins(), |
| 629 | fVtxAxis.GetXmin(), |
| 630 | fVtxAxis.GetXmax()); |
| 563a673f |
631 | selEff->Sumw2(); |
| 632 | selEff->SetDirectory(0); |
| 633 | selEff->SetFillColor(kMagenta+1); |
| 634 | selEff->SetFillStyle(3001); |
| 635 | selEff->Add(eventsAll); |
| 636 | selEff->Divide(eventsTrVtx); |
| 637 | list->Add(selEff); |
| 638 | |
| 639 | // Loop over vertex bins and do vertex dependendt stuff |
| 640 | for (Int_t v = 1; v <= fVtxAxis.GetNbins(); v++) { |
| 641 | // Make a sub-list in the output |
| 642 | TList* vl = new TList; |
| 0bd4b00f |
643 | vl->SetName(Form("vtx%02d", v)); |
| 563a673f |
644 | list->Add(vl); |
| 645 | |
| 646 | // Get event counts |
| 0bd4b00f |
647 | Int_t nEventsAll = eventsAll->GetBinContent(v); |
| 648 | Int_t nEventsTr = eventsTr->GetBinContent(v); |
| 649 | Int_t nEventsVtx = eventsVtx->GetBinContent(v); |
| 650 | Int_t nEventsTrVtx = eventsTrVtx->GetBinContent(v); |
| 563a673f |
651 | |
| 652 | // Project event histograms, set names, and store |
| 653 | TH2D* primIAllV = GetVertexProj(v, primIAll); |
| 654 | TH2D* primOAllV = GetVertexProj(v, primOAll); |
| 655 | TH2D* primITrVtxV = GetVertexProj(v, primITrVtx); |
| 656 | TH2D* primOTrVtxV = GetVertexProj(v, primOTrVtx); |
| 657 | vl->Add(primIAllV); |
| 658 | vl->Add(primOAllV); |
| 659 | vl->Add(primITrVtxV); |
| 660 | vl->Add(primOTrVtxV); |
| 661 | |
| 0bd4b00f |
662 | primIAllV->Scale(1. / nEventsAll); |
| 663 | primOAllV->Scale(1. / nEventsAll); |
| 664 | primITrVtxV->Scale(1. / nEventsTr); |
| 665 | primOTrVtxV->Scale(1. / nEventsTr); |
| 563a673f |
666 | |
| 667 | // Calculate the vertex bias on the d^2N/detadphi |
| 0bd4b00f |
668 | TH2D* selBiasI = |
| 669 | static_cast<TH2D*>(primITrVtxV->Clone(Form("selBiasI%02d",v))); |
| 670 | TH2D* selBiasO = |
| 671 | static_cast<TH2D*>(primOTrVtxV->Clone(Form("selBiasO%02d",v))); |
| 563a673f |
672 | selBiasI->SetTitle(Form("Event selection bias in vertex bin %d", v)); |
| 673 | selBiasO->SetTitle(Form("Event selection bias in vertex bin %d", v)); |
| 674 | selBiasI->Divide(primIAllV); |
| 675 | selBiasO->Divide(primOAllV); |
| 676 | vl->Add(selBiasI); |
| 677 | vl->Add(selBiasO); |
| 678 | |
| 679 | for (UShort_t d = 1; d <= 3; d++) { |
| 680 | UShort_t nQ = (d == 1 ? 1 : 2); |
| 681 | for (UShort_t q = 0; q < nQ; q++) { |
| 682 | Char_t r = (q == 0 ? 'I' : 'O'); |
| 683 | TH3D* hits3D = 0; |
| 684 | TH2D* prim2D = (q == 0 ? primITrVtxV : primOTrVtxV); |
| 685 | switch (d) { |
| 686 | case 1: hits3D = hits1i; break; |
| 687 | case 2: hits3D = (q == 0 ? hits2i : hits2o); break; |
| 688 | case 3: hits3D = (q == 0 ? hits3i : hits3o); break; |
| 689 | } |
| 690 | |
| 691 | TH2D* sec = GetVertexProj(v, hits3D); |
| 692 | sec->SetName(Form("secondaryFMD%d%c_vtx%02d", d, r, v)); |
| 693 | sec->SetTitle(Form("Secondary correction map for FMD%d%c " |
| 694 | "in vertex bin %d", d, r, v)); |
| 695 | sec->Divide(prim2D); |
| 696 | vl->Add(sec); |
| 697 | |
| 698 | if (v > 1) continue; |
| 699 | |
| 700 | // Do the double hit correction (only done once per ring in |
| 701 | // the vertex loop) |
| 0bd4b00f |
702 | TH1D* hStrips = 0; |
| 563a673f |
703 | switch (d) { |
| 0bd4b00f |
704 | case 1: hStrips = strips1i; break; |
| 705 | case 2: hStrips = (q == 0 ? strips2i : strips2o); break; |
| 706 | case 3: hStrips = (q == 0 ? strips3i : strips3o); break; |
| 563a673f |
707 | } |
| 708 | |
| 709 | TH2D* hits2D = GetVertexProj(v, hits3D); |
| 710 | TH1D* doubleHit = hits2D->ProjectionX(Form("doubleHitFMD%d%c",d,r)); |
| 711 | doubleHit->SetTitle(Form("Double hit correction for FMD%d%c",d,r)); |
| 712 | doubleHit->SetDirectory(0); |
| 713 | doubleHit->SetFillColor(kGreen+1); |
| 714 | doubleHit->SetFillStyle(3001); |
| 715 | doubleHit->Sumw2(); |
| 0bd4b00f |
716 | doubleHit->Divide(hStrips); |
| 563a673f |
717 | list->Add(doubleHit); |
| 718 | } |
| 719 | } |
| 720 | } |
| 721 | } |
| 722 | |
| 723 | //____________________________________________________________________ |
| 724 | void |
| 0bd4b00f |
725 | AliForwardMCCorrectionsTask::Print(Option_t*) const |
| 563a673f |
726 | { |
| 727 | } |
| 728 | |
| 729 | // |
| 730 | // EOF |
| 731 | // |
git://git.uio.no / u / mrichter / AliRoot.git / blame