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