]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliBasedNdetaTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Various fixes to deal with centrality
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliBasedNdetaTask.cxx
CommitLineData
ce85db45 1//====================================================================
2#include "AliBasedNdetaTask.h"
3#include <TMath.h>
4#include <TH2D.h>
5#include <TH1D.h>
6#include <THStack.h>
7#include <TList.h>
e58000b7 8#include <TParameter.h>
ce85db45 9#include <AliAnalysisManager.h>
10#include <AliAODEvent.h>
11#include <AliAODHandler.h>
12#include <AliAODInputHandler.h>
13#include "AliForwardUtil.h"
14#include "AliAODForwardMult.h"
3846ec25 15#include <TFile.h>
ce85db45 16
17//____________________________________________________________________
18AliBasedNdetaTask::AliBasedNdetaTask()
19 : AliAnalysisTaskSE(),
ce85db45 20 fSums(0), // Container of sums
e1f47419 21 fOutput(0), // Container of output
ce85db45 22 fVtxMin(0), // Minimum v_z
23 fVtxMax(0), // Maximum v_z
24 fTriggerMask(0), // Trigger mask
25 fRebin(0), // Rebinning factor
26 fCutEdges(false),
27 fSymmetrice(true),
e58000b7 28 fCorrEmpty(true),
29 fTriggerEff(1),
3846ec25 30 fShapeCorr(0),
e1f47419 31 fListOfCentralities(0),
32 fUseShapeCorr(true),
33 fSNNString(0),
34 fSysString(0),
e308a636 35 fCent(0),
36 fCentAxis(0)
fb3430ac 37{
38 //
39 // Constructor
40 //
41}
ce85db45 42
43//____________________________________________________________________
44AliBasedNdetaTask::AliBasedNdetaTask(const char* name)
45 : AliAnalysisTaskSE(name),
ce85db45 46 fSums(0), // Container of sums
e1f47419 47 fOutput(0), // Container of output
ce85db45 48 fVtxMin(-10), // Minimum v_z
49 fVtxMax(10), // Maximum v_z
50 fTriggerMask(AliAODForwardMult::kInel),
51 fRebin(5), // Rebinning factor
52 fCutEdges(false),
53 fSymmetrice(true),
e58000b7 54 fCorrEmpty(true),
55 fTriggerEff(1),
3846ec25 56 fShapeCorr(0),
e1f47419 57 fListOfCentralities(0),
58 fUseShapeCorr(true),
59 fSNNString(0),
60 fSysString(0),
e308a636 61 fCent(0),
62 fCentAxis(0)
ce85db45 63{
fb3430ac 64 //
65 // Constructor
66 //
e308a636 67 fListOfCentralities = new TObjArray(1);
fb3430ac 68
ce85db45 69 // Output slot #1 writes into a TH1 container
70 DefineOutput(1, TList::Class());
71 DefineOutput(2, TList::Class());
72}
73
74//____________________________________________________________________
75AliBasedNdetaTask::AliBasedNdetaTask(const AliBasedNdetaTask& o)
76 : AliAnalysisTaskSE(o),
ce85db45 77 fSums(o.fSums), // TList* - Container of sums
e1f47419 78 fOutput(o.fOutput), // Container of output
ce85db45 79 fVtxMin(o.fVtxMin), // Double_t - Minimum v_z
80 fVtxMax(o.fVtxMax), // Double_t - Maximum v_z
81 fTriggerMask(o.fTriggerMask),// Int_t - Trigger mask
82 fRebin(o.fRebin), // Int_t - Rebinning factor
83 fCutEdges(o.fCutEdges), // Bool_t - Whether to cut edges when rebinning
84 fSymmetrice(true),
e58000b7 85 fCorrEmpty(true),
86 fTriggerEff(o.fTriggerEff),
3846ec25 87 fShapeCorr(o.fShapeCorr),
e1f47419 88 fListOfCentralities(o.fListOfCentralities),
89 fUseShapeCorr(o.fUseShapeCorr),
90 fSNNString(o.fSNNString),
91 fSysString(o.fSysString),
e308a636 92 fCent(o.fCent),
93 fCentAxis(o.fCentAxis)
ce85db45 94{}
95
96//____________________________________________________________________
97AliBasedNdetaTask::~AliBasedNdetaTask()
98{
fb3430ac 99 //
100 // Destructor
101 //
ce85db45 102 if (fSums) {
103 fSums->Delete();
104 delete fSums;
105 fSums = 0;
106 }
e1f47419 107 if (fOutput) {
108 fOutput->Delete();
109 delete fOutput;
110 fOutput = 0;
ce85db45 111 }
112}
113
e1f47419 114//________________________________________________________________________
115void
e308a636 116AliBasedNdetaTask::SetCentralityAxis(UShort_t n, Short_t* bins)
117{
118 if (!fCentAxis) {
119 fCentAxis = new TAxis();
120 fCentAxis->SetName("centAxis");
121 fCentAxis->SetTitle("Centrality [%]");
122 }
123 TArrayD dbins(n+1);
124 for (UShort_t i = 0; i <= n; i++)
125 dbins[i] = (bins[i] == 100 ? 100.1 : bins[i]);
126 fCentAxis->Set(n, dbins.GetArray());
127}
128
129//________________________________________________________________________
130void
131AliBasedNdetaTask::AddCentralityBin(UShort_t at, Short_t low, Short_t high)
e1f47419 132{
133 //
134 // Add a centrality bin
135 //
136 // Parameters:
137 // low Low cut
138 // high High cut
139 //
e308a636 140 CentralityBin* bin = MakeCentralityBin(GetName(), low, high);
141 AliInfo(Form("Adding centrality bin %p [%3d,%3d] @ %d", bin, low, high, at));
142 fListOfCentralities->AddAtAndExpand(bin, at);
e1f47419 143}
144
145//________________________________________________________________________
146AliBasedNdetaTask::CentralityBin*
147AliBasedNdetaTask::MakeCentralityBin(const char* name,
148 Short_t low, Short_t high) const
149{
150 //
151 // Make a centrality bin
152 //
153 // Parameters:
154 // name Name used for histograms
155 // low Low cut in percent
156 // high High cut in percent
157 //
158 // Return:
159 // A newly created centrality bin
160 //
161 return new CentralityBin(name, low, high);
162}
ce85db45 163//________________________________________________________________________
164void
165AliBasedNdetaTask::SetTriggerMask(const char* mask)
166{
fb3430ac 167 //
168 // Set the trigger maskl
169 //
170 // Parameters:
171 // mask Trigger mask
172 //
ce85db45 173 UShort_t trgMask = 0;
174 TString trgs(mask);
175 trgs.ToUpper();
176 TObjString* trg;
177 TIter next(trgs.Tokenize(" ,|"));
178 while ((trg = static_cast<TObjString*>(next()))) {
179 TString s(trg->GetString());
180 if (s.IsNull()) continue;
181 if (s.CompareTo("INEL") == 0) trgMask = AliAODForwardMult::kInel;
182 else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
183 else if (s.CompareTo("NSD") == 0) trgMask = AliAODForwardMult::kNSD;
184 else
185 Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
186 }
187 if (trgMask == 0) trgMask = 1;
188 SetTriggerMask(trgMask);
189}
190
e58000b7 191//________________________________________________________________________
192void
193AliBasedNdetaTask::SetShapeCorrection(const TH1* c)
194{
fb3430ac 195 //
196 // Set the shape correction (a.k.a., track correction) for selected
197 // trigger(s)
198 //
199 // Parameters:
200 // h Correction
201 //
e58000b7 202 if (!c) return;
203 fShapeCorr = static_cast<TH1*>(c->Clone());
204 fShapeCorr->SetDirectory(0);
205}
206
ce85db45 207//________________________________________________________________________
208void
209AliBasedNdetaTask::UserCreateOutputObjects()
210{
fb3430ac 211 //
212 // Create output objects.
213 //
214 // This is called once per slave process
215 //
ce85db45 216 fSums = new TList;
217 fSums->SetName(Form("%s_sums", GetName()));
218 fSums->SetOwner();
219
e1f47419 220 // Automatically add 'all' centrality bin if nothing has been defined.
e308a636 221 AddCentralityBin(0, 0, 0);
222 if (fCentAxis && fCentAxis->GetNbins() > 0 && fCentAxis->GetXbins()) {
223 const TArrayD* bins = fCentAxis->GetXbins();
224 Int_t nbin = fCentAxis->GetNbins();
225 for (Int_t i = 0; i < nbin; i++)
226 AddCentralityBin(i+1, Short_t((*bins)[i]), Short_t((*bins)[i+1]));
227 }
228 fSums->Add(fCentAxis);
229
ce85db45 230
e1f47419 231 // Centrality histogram
232 fCent = new TH1D("cent", "Centrality", 100, 0, 100);
233 fCent->SetDirectory(0);
234 fCent->SetXTitle(0);
235 fSums->Add(fCent);
236
237 // Loop over centrality bins
238 TIter next(fListOfCentralities);
239 CentralityBin* bin = 0;
240 while ((bin = static_cast<CentralityBin*>(next())))
241 bin->CreateOutputObjects(fSums);
ce85db45 242
243 // Check that we have an AOD input handler
244 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
245 AliAODInputHandler* ah =
246 dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
247 if (!ah) AliFatal("No AOD input handler set in analysis manager");
248
249 // Post data for ALL output slots >0 here, to get at least an empty histogram
250 PostData(1, fSums);
251}
ce85db45 252//____________________________________________________________________
253void
254AliBasedNdetaTask::UserExec(Option_t *)
255{
fb3430ac 256 //
257 // Process a single event
258 //
259 // Parameters:
260 // option Not used
261 //
ce85db45 262 // Main loop
263 AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
264 if (!aod) {
265 AliError("Cannot get the AOD event");
266 return;
267 }
e58000b7 268
e1f47419 269 TObject* obj = aod->FindListObject("Forward");
270 if (!obj) {
271 AliWarning("No forward object found");
272 return;
678a4979 273 }
e1f47419 274 AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
e308a636 275
e1f47419 276 // Fill centrality histogram
e308a636 277 Float_t cent = forward->GetCentrality();
278 fCent->Fill(cent);
e1f47419 279
ce85db45 280 // Get the histogram(s)
e58000b7 281 TH2D* data = GetHistogram(aod, false);
ce85db45 282 TH2D* dataMC = GetHistogram(aod, true);
283
e1f47419 284 // Loop over centrality bins
e308a636 285 CentralityBin* allBin = static_cast<CentralityBin*>(fListOfCentralities->At(0));
286 allBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
287
288 // Find this centrality bin
289 if (fCentAxis && fCentAxis->GetNbins() > 0) {
290 Int_t icent = fCentAxis->FindBin(cent);
291 CentralityBin* thisBin = 0;
292 if (icent >= 1 && icent <= fCentAxis->GetNbins())
293 thisBin = static_cast<CentralityBin*>(fListOfCentralities->At(icent));
294 if (thisBin)
295 thisBin->ProcessEvent(forward, fTriggerMask, fVtxMin, fVtxMax, data, dataMC);
296 }
ce85db45 297
e1f47419 298 // Here, we get the update
299 if (!fSNNString) {
300 UShort_t sNN = forward->GetSNN();
301 fSNNString = new TNamed("sNN", "");
302 fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
303 fSNNString->SetUniqueID(sNN);
304 fSums->Add(fSNNString);
e58000b7 305
e1f47419 306 UShort_t sys = forward->GetSystem();
307 fSysString = new TNamed("sys", "");
308 fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
309 fSysString->SetUniqueID(sys);
310 fSums->Add(fSysString);
3846ec25 311 }
e58000b7 312
ce85db45 313 PostData(1, fSums);
314}
315
316//________________________________________________________________________
317void
318AliBasedNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
319 const char* title, const char* ytitle)
320{
fb3430ac 321 //
322 // Set histogram graphical options, etc.
323 //
324 // Parameters:
325 // h Histogram to modify
326 // colour Marker color
327 // marker Marker style
328 // title Title of histogram
329 // ytitle Title on y-axis.
330 //
ce85db45 331 h->SetTitle(title);
332 h->SetMarkerColor(colour);
333 h->SetMarkerStyle(marker);
334 h->SetMarkerSize(1);
335 h->SetFillStyle(0);
336 h->SetYTitle(ytitle);
337 h->GetXaxis()->SetTitleFont(132);
338 h->GetXaxis()->SetLabelFont(132);
339 h->GetXaxis()->SetNdivisions(10);
340 h->GetYaxis()->SetTitleFont(132);
341 h->GetYaxis()->SetLabelFont(132);
342 h->GetYaxis()->SetNdivisions(10);
343 h->GetYaxis()->SetDecimals();
344 h->SetStats(0);
345}
346
347//________________________________________________________________________
348TH1D*
349AliBasedNdetaTask::ProjectX(const TH2D* h,
350 const char* name,
351 Int_t firstbin,
352 Int_t lastbin,
353 bool corr,
e1f47419 354 bool error)
ce85db45 355{
fb3430ac 356 //
357 // Project onto the X axis
358 //
359 // Parameters:
360 // h 2D histogram
361 // name New name
362 // firstbin First bin to use
363 // lastbin Last bin to use
364 // error Whether to calculate errors
365 //
366 // Return:
367 // Newly created histogram or null
368 //
ce85db45 369 if (!h) return 0;
e1f47419 370#if USE_ROOT_PROJECT
371 return h->ProjectionX(name, firstbin, lastbin, (error ? "e" : ""));
372#endif
ce85db45 373
374 TAxis* xaxis = h->GetXaxis();
375 TAxis* yaxis = h->GetYaxis();
376 TH1D* ret = new TH1D(name, h->GetTitle(), xaxis->GetNbins(),
377 xaxis->GetXmin(), xaxis->GetXmax());
378 static_cast<const TAttLine*>(h)->Copy(*ret);
379 static_cast<const TAttFill*>(h)->Copy(*ret);
380 static_cast<const TAttMarker*>(h)->Copy(*ret);
381 ret->GetXaxis()->ImportAttributes(xaxis);
382
383 Int_t first = firstbin;
384 Int_t last = lastbin;
385 if (first < 0) first = 0;
386 else if (first >= yaxis->GetNbins()+1) first = yaxis->GetNbins();
387 if (last < 0) last = yaxis->GetNbins();
388 else if (last > yaxis->GetNbins()+1) last = yaxis->GetNbins();
389 if (last-first < 0) {
e1f47419 390 AliWarningGeneral("AliBasedNdetaTask",
391 Form("Nothing to project [%d,%d]", first, last));
ce85db45 392 return 0;
393
394 }
395
396 // Loop over X bins
397 // AliInfo(Form("Projecting bins [%d,%d] of %s", first, last, h->GetName()));
398 Int_t ybins = (last-first+1);
399 for (Int_t xbin = 0; xbin <= xaxis->GetNbins()+1; xbin++) {
400 Double_t content = 0;
401 Double_t error2 = 0;
402 Int_t nbins = 0;
403
404
405 for (Int_t ybin = first; ybin <= last; ybin++) {
406 Double_t c1 = h->GetCellContent(xbin, ybin);
407 Double_t e1 = h->GetCellError(xbin, ybin);
408
409 // Ignore empty bins
410 if (c1 < 1e-12) continue;
411 if (e1 < 1e-12) {
412 if (error) continue;
413 e1 = 1;
414 }
415
416 content += c1;
417 error2 += e1*e1;
418 nbins++;
419 } // for (ybin)
420 if(content > 0 && nbins > 0) {
421 Double_t factor = (corr ? Double_t(ybins) / nbins : 1);
422 if (error) {
423 // Calculate weighted average
424 ret->SetBinContent(xbin, content * factor);
425 ret->SetBinError(xbin, factor * TMath::Sqrt(error2));
426 }
427 else
428 ret->SetBinContent(xbin, factor * content);
429 }
430 } // for (xbin)
431
432 return ret;
433}
434
435//________________________________________________________________________
436void
437AliBasedNdetaTask::Terminate(Option_t *)
438{
fb3430ac 439 //
440 // Called at end of event processing..
441 //
442 // This is called once in the master
443 //
444 // Parameters:
445 // option Not used
446 //
ce85db45 447 // Draw result to screen, or perform fitting, normalizations Called
448 // once at the end of the query
449
450 fSums = dynamic_cast<TList*> (GetOutputData(1));
451 if(!fSums) {
452 AliError("Could not retrieve TList fSums");
453 return;
454 }
455
e1f47419 456 fOutput = new TList;
457 fOutput->SetName(Form("%s_result", GetName()));
458 fOutput->SetOwner();
e58000b7 459
e1f47419 460 fSNNString = static_cast<TNamed*>(fSums->FindObject("sNN"));
461 fSysString = static_cast<TNamed*>(fSums->FindObject("sys"));
e308a636 462 fCentAxis = static_cast<TAxis*>(fSums->FindObject("centAxis"));
ce85db45 463
e1f47419 464 if(fSysString && fSNNString &&
465 fSysString->GetUniqueID() == AliForwardUtil::kPP)
466 LoadNormalizationData(fSysString->GetUniqueID(),
467 fSNNString->GetUniqueID());
ce85db45 468
e1f47419 469 // Loop over centrality bins
470 TIter next(fListOfCentralities);
471 CentralityBin* bin = 0;
472 while ((bin = static_cast<CentralityBin*>(next())))
473 bin->End(fSums, fOutput, fUseShapeCorr ? fShapeCorr : 0, fTriggerEff,
474 fSymmetrice, fRebin, fCorrEmpty, fCutEdges,
475 fVtxMin, fVtxMax, fTriggerMask);
ce85db45 476
e1f47419 477 // Output collision energy string
478 if (fSNNString) fOutput->Add(fSNNString->Clone());
ce85db45 479
e1f47419 480 // Output collision system string
481 if (fSysString) fOutput->Add(fSysString->Clone());
ce85db45 482
e308a636 483 // Output centrality axis
484 if (fCentAxis) fOutput->Add(fCentAxis);
485
e1f47419 486 // Output trigger string
ce85db45 487 TNamed* trigString =
488 new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
489 trigString->SetUniqueID(fTriggerMask);
490 fOutput->Add(trigString);
491
e1f47419 492 // Output vertex axis
ce85db45 493 TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
494 vtxAxis->SetName("vtxAxis");
495 vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
496 fOutput->Add(vtxAxis);
e1f47419 497
498 // Output trigger efficiency and shape correction
e58000b7 499 fOutput->Add(new TParameter<Double_t>("triggerEff", fTriggerEff));
500 if (fShapeCorr) fOutput->Add(fShapeCorr);
501
ce85db45 502 PostData(2, fOutput);
503}
3846ec25 504//________________________________________________________________________
505void
506AliBasedNdetaTask::LoadNormalizationData(UShort_t sys, UShort_t energy)
507{
e1f47419 508 // Load the normalisation data for dN/deta for pp INEL, INEL>0, and NSD
3846ec25 509 TString type("pp");
3846ec25 510 TString snn("900");
511 if(energy == 7000) snn.Form("7000");
512 if(energy == 2750) snn.Form("2750");
513
e1f47419 514 if(fShapeCorr && (fTriggerEff != 1)) {
515 AliInfo("Objects already set for normalization - no action taken");
516 return;
517 }
3846ec25 518
e1f47419 519 TFile* fin = TFile::Open(Form("$ALICE_ROOT/PWG2/FORWARD/corrections/"
520 "Normalization/normalizationHists_%s_%s.root",
521 type.Data(),snn.Data()));
522 if(!fin) {
523 AliWarning(Form("no file for normalization of %d/%d", sys, energy));
524 return;
525 }
3846ec25 526
e1f47419 527 // Shape correction
528 TString shapeCorName(Form("h%sNormalization",
529 fTriggerMask == AliAODForwardMult::kInel ? "Inel" :
530 fTriggerMask == AliAODForwardMult::kNSD ? "NSD" :
531 fTriggerMask == AliAODForwardMult::kInelGt0 ?
532 "InelGt0" : "All"));
533 TH2F* shapeCor = dynamic_cast<TH2F*>(fin->Get(shapeCorName));
534 if (shapeCor) SetShapeCorrection(shapeCor);
ce85db45 535
e1f47419 536 // Trigger efficiency
537 TString effName(Form("%sTriggerEff",
538 fTriggerMask == AliAODForwardMult::kInel ? "inel" :
539 fTriggerMask == AliAODForwardMult::kNSD ? "nsd" :
540 fTriggerMask == AliAODForwardMult::kInelGt0 ?
541 "inelgt0" : "all"));
542 TParameter<float>* eff = static_cast<TParameter<float>*>(fin->Get(effName));
543 if (eff) SetTriggerEff(eff->GetVal());
3846ec25 544
e1f47419 545 // Print - out
546 if (shapeCor && eff) AliInfo("Loaded objects for normalization.");
3846ec25 547}
ce85db45 548//________________________________________________________________________
549TH1D*
e1f47419 550AliBasedNdetaTask::Rebin(const TH1D* h, Int_t rebin, Bool_t cutEdges)
ce85db45 551{
fb3430ac 552 //
553 // Make a copy of the input histogram and rebin that histogram
554 //
555 // Parameters:
556 // h Histogram to rebin
557 //
558 // Return:
559 // New (rebinned) histogram
560 //
e1f47419 561 if (rebin <= 1) return 0;
ce85db45 562
563 Int_t nBins = h->GetNbinsX();
e1f47419 564 if(nBins % rebin != 0) {
565 AliWarningGeneral("AliBasedNdetaTask",
566 Form("Rebin factor %d is not a devisor of current number "
567 "of bins %d in the histogram %s",
568 rebin, nBins, h->GetName()));
ce85db45 569 return 0;
570 }
571
572 // Make a copy
573 TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
e1f47419 574 h->GetName(), rebin)));
575 tmp->Rebin(rebin);
ce85db45 576 tmp->SetDirectory(0);
577
578 // The new number of bins
e1f47419 579 Int_t nBinsNew = nBins / rebin;
ce85db45 580 for(Int_t i = 1;i<= nBinsNew; i++) {
581 Double_t content = 0;
582 Double_t sumw = 0;
583 Double_t wsum = 0;
584 Int_t nbins = 0;
e1f47419 585 for(Int_t j = 1; j<=rebin;j++) {
586 Int_t bin = (i-1)*rebin + j;
ce85db45 587 Double_t c = h->GetBinContent(bin);
ce85db45 588 if (c <= 0) continue;
678a4979 589
e1f47419 590 if (cutEdges) {
ce85db45 591 if (h->GetBinContent(bin+1)<=0 ||
3846ec25 592 h->GetBinContent(bin-1)<=0) {
e1f47419 593#if 0
594 AliWarningGeneral("AliBasedNdetaTask",
595 Form("removing bin %d=%f of %s (%d=%f,%d=%f)",
596 bin, c, h->GetName(),
597 bin+1, h->GetBinContent(bin+1),
598 bin-1, h->GetBinContent(bin-1)));
599#endif
ce85db45 600 continue;
601 }
602 }
603 Double_t e = h->GetBinError(bin);
604 Double_t w = 1 / (e*e); // 1/c/c
605 content += c;
606 sumw += w;
607 wsum += w * c;
608 nbins++;
609 }
610
611 if(content > 0 && nbins > 0) {
612 tmp->SetBinContent(i, wsum / sumw);
613 tmp->SetBinError(i,1./TMath::Sqrt(sumw));
614 }
615 }
616
617 return tmp;
618}
619
620//__________________________________________________________________
ce85db45 621TH1*
e1f47419 622AliBasedNdetaTask::Symmetrice(const TH1* h)
ce85db45 623{
fb3430ac 624 //
625 // Make an extension of @a h to make it symmetric about 0
626 //
627 // Parameters:
628 // h Histogram to symmertrice
629 //
630 // Return:
631 // Symmetric extension of @a h
632 //
ce85db45 633 Int_t nBins = h->GetNbinsX();
634 TH1* s = static_cast<TH1*>(h->Clone(Form("%s_mirror", h->GetName())));
635 s->SetTitle(Form("%s (mirrored)", h->GetTitle()));
636 s->Reset();
637 s->SetBins(nBins, -h->GetXaxis()->GetXmax(), -h->GetXaxis()->GetXmin());
638 s->SetMarkerColor(h->GetMarkerColor());
639 s->SetMarkerSize(h->GetMarkerSize());
640 s->SetMarkerStyle(h->GetMarkerStyle()+4);
641 s->SetFillColor(h->GetFillColor());
642 s->SetFillStyle(h->GetFillStyle());
643 s->SetDirectory(0);
644
645 // Find the first and last bin with data
646 Int_t first = nBins+1;
647 Int_t last = 0;
648 for (Int_t i = 1; i <= nBins; i++) {
649 if (h->GetBinContent(i) <= 0) continue;
650 first = TMath::Min(first, i);
651 last = TMath::Max(last, i);
652 }
653
654 Double_t xfirst = h->GetBinCenter(first-1);
655 Int_t f1 = h->GetXaxis()->FindBin(-xfirst);
656 Int_t l2 = s->GetXaxis()->FindBin(xfirst);
657 for (Int_t i = f1, j=l2; i <= last; i++,j--) {
658 s->SetBinContent(j, h->GetBinContent(i));
659 s->SetBinError(j, h->GetBinError(i));
660 }
661 // Fill in overlap bin
662 s->SetBinContent(l2+1, h->GetBinContent(first));
663 s->SetBinError(l2+1, h->GetBinError(first));
664 return s;
665}
e1f47419 666
667//====================================================================
668AliBasedNdetaTask::CentralityBin::CentralityBin()
669 : TNamed("", ""),
670 fSums(0),
671 fOutput(0),
672 fSum(0),
673 fSumMC(0),
674 fTriggers(0),
675 fLow(0),
676 fHigh(0)
677{
678 //
679 // Constructor
680 //
681}
682//____________________________________________________________________
683AliBasedNdetaTask::CentralityBin::CentralityBin(const char* name,
684 Short_t low, Short_t high)
685 : TNamed(name, ""),
686 fSums(0),
687 fOutput(0),
688 fSum(0),
689 fSumMC(0),
690 fTriggers(0),
691 fLow(low),
692 fHigh(high)
693{
694 //
695 // Constructor
696 //
697 // Parameters:
698 // name Name used for histograms (e.g., Forward)
699 // low Lower centrality cut in percent
700 // high Upper centrality cut in percent
701 //
702 if (low <= 0 && high <= 0) {
703 fLow = 0;
704 fHigh = 0;
705 SetTitle("All centralities");
706 }
707 else {
708 fLow = low;
709 fHigh = high;
710 SetTitle(Form("Centrality bin from %3d%% to %3d%%", low, high));
711 }
712}
713//____________________________________________________________________
714AliBasedNdetaTask::CentralityBin::CentralityBin(const CentralityBin& o)
715 : TNamed(o),
716 fSums(o.fSums),
717 fOutput(o.fOutput),
718 fSum(o.fSum),
719 fSumMC(o.fSumMC),
720 fTriggers(o.fTriggers),
721 fLow(o.fLow),
722 fHigh(o.fHigh)
723{
724 //
725 // Copy constructor
726 //
727 // Parameters:
728 // other Object to copy from
729 //
730}
731//____________________________________________________________________
732AliBasedNdetaTask::CentralityBin::~CentralityBin()
733{
734 //
735 // Destructor
736 //
737 if (fSums) fSums->Delete();
738 if (fOutput) fOutput->Delete();
739}
740
741//____________________________________________________________________
742AliBasedNdetaTask::CentralityBin&
743AliBasedNdetaTask::CentralityBin::operator=(const CentralityBin& o)
744{
745 //
746 // Assignment operator
747 //
748 // Parameters:
749 // other Object to assign from
750 //
751 // Return:
752 // Reference to this
753 //
754 SetName(o.GetName());
755 SetTitle(o.GetTitle());
756 fSums = o.fSums;
757 fOutput = o.fOutput;
758 fSum = o.fSum;
759 fSumMC = o.fSumMC;
760 fTriggers = o.fTriggers;
761 fLow = o.fLow;
762 fHigh = o.fHigh;
763
764 return *this;
765}
766//____________________________________________________________________
767const char*
768AliBasedNdetaTask::CentralityBin::GetListName() const
769{
770 //
771 // Get the list name
772 //
773 // Return:
774 // List Name
775 //
776 if (IsAllBin()) return "all";
777 return Form("cent%03d_%03d", fLow, fHigh);
778}
779//____________________________________________________________________
780void
781AliBasedNdetaTask::CentralityBin::CreateOutputObjects(TList* dir)
782{
783 //
784 // Create output objects
785 //
786 // Parameters:
787 // dir Parent list
788 //
789 fSums = new TList;
790 fSums->SetName(GetListName());
791 fSums->SetOwner();
792 dir->Add(fSums);
793
794 fTriggers = new TH1D("triggers", "Number of triggers",
795 kAccepted, 1, kAccepted);
796 fTriggers->SetYTitle("# of events");
797 fTriggers->GetXaxis()->SetBinLabel(kAll, "All events");
798 fTriggers->GetXaxis()->SetBinLabel(kB, "w/B trigger");
799 fTriggers->GetXaxis()->SetBinLabel(kA, "w/A trigger");
800 fTriggers->GetXaxis()->SetBinLabel(kC, "w/C trigger");
801 fTriggers->GetXaxis()->SetBinLabel(kE, "w/E trigger");
802 fTriggers->GetXaxis()->SetBinLabel(kMB, "w/Collision trigger");
803 fTriggers->GetXaxis()->SetBinLabel(kWithVertex, "w/Vertex");
804 fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
805 fTriggers->GetXaxis()->SetBinLabel(kPileUp, "w/Pileup");
806 fTriggers->GetXaxis()->SetBinLabel(kAccepted, "Accepted by cut");
807 fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
808 fTriggers->SetFillColor(kRed+1);
809 fTriggers->SetFillStyle(3001);
810 fTriggers->SetStats(0);
811 fSums->Add(fTriggers);
812}
813//____________________________________________________________________
814void
815AliBasedNdetaTask::CentralityBin::CreateSums(const TH2D* data, const TH2D* mc)
816{
817 //
818 // Create sum histogram
819 //
820 // Parameters:
821 // data Data histogram to clone
822 // mc (optional) MC histogram to clone
823 //
824 fSum = static_cast<TH2D*>(data->Clone(GetName()));
825 fSum->SetDirectory(0);
826 fSum->Reset();
827 fSums->Add(fSum);
828
829 // If no MC data is given, then do not create MC sum histogram
830 if (!mc) return;
831
832 fSumMC = static_cast<TH2D*>(mc->Clone(Form("%sMC", GetName())));
833 fSumMC->SetDirectory(0);
834 fSumMC->Reset();
835 fSums->Add(fSumMC);
836}
837
838//____________________________________________________________________
839Bool_t
840AliBasedNdetaTask::CentralityBin::CheckEvent(const AliAODForwardMult* forward,
841 Int_t triggerMask,
842 Double_t vzMin, Double_t vzMax)
843{
844 //
845 // Check the trigger, vertex, and centrality
846 //
847 // Parameters:
848 // aod Event input
849 //
850 // Return:
851 // true if the event is to be used
852 //
853 if (!forward) return false;
854
855 // Check the centrality class unless this is the 'all' bin
e308a636 856 // if (!IsAllBin()) {
857 // Double_t centrality = forward->GetCentrality();
858 // if (centrality < fLow || centrality >= fHigh) return false;
859 // }
e1f47419 860
861 fTriggers->AddBinContent(kAll);
862 if (forward->IsTriggerBits(AliAODForwardMult::kB))
863 fTriggers->AddBinContent(kB);
864 if (forward->IsTriggerBits(AliAODForwardMult::kA))
865 fTriggers->AddBinContent(kA);
866 if (forward->IsTriggerBits(AliAODForwardMult::kC))
867 fTriggers->AddBinContent(kC);
868 if (forward->IsTriggerBits(AliAODForwardMult::kE))
869 fTriggers->AddBinContent(kE);
870 if (forward->IsTriggerBits(AliAODForwardMult::kInel))
871 fTriggers->AddBinContent(kMB);
872 if (forward->IsTriggerBits(AliAODForwardMult::kMCNSD))
873 fTriggers->AddBinContent(kMCNSD);
874 if (forward->IsTriggerBits(AliAODForwardMult::kPileUp))
875 fTriggers->AddBinContent(kPileUp);
876
877 // Check if we have an event of interest.
878 if (!forward->IsTriggerBits(triggerMask)) return false;
879
880 //Check for pileup
881 if (forward->IsTriggerBits(AliAODForwardMult::kPileUp)) return false;
882
883 fTriggers->AddBinContent(kWithTrigger);
884
885 // Check that we have a valid vertex
886 if (!forward->HasIpZ()) return false;
887 fTriggers->AddBinContent(kWithVertex);
888
889 // Check that vertex is within cuts
890 if (!forward->InRange(vzMin, vzMax)) return false;
891 fTriggers->AddBinContent(kAccepted);
892
893 return true;
894}
895
896
897//____________________________________________________________________
898void
899AliBasedNdetaTask::CentralityBin::ProcessEvent(const AliAODForwardMult* forward,
900 Int_t triggerMask,
901 Double_t vzMin, Double_t vzMax,
902 const TH2D* data, const TH2D* mc)
903{
904 //
905 // Process an event
906 //
907 // Parameters:
908 // forward Forward data (for trigger, vertex, & centrality)
909 // triggerMask Trigger mask
910 // vzMin Minimum IP z coordinate
911 // vzMax Maximum IP z coordinate
912 // data Data histogram
913 // mc MC histogram
914 //
915 if (!CheckEvent(forward, triggerMask, vzMin, vzMax)) return;
916 if (!data) return;
917 if (!fSum) CreateSums(data, mc);
918
919 fSum->Add(data);
920 if (mc) fSumMC->Add(mc);
921}
922
923//________________________________________________________________________
924void
925AliBasedNdetaTask::CentralityBin::End(TList* sums,
926 TList* results,
927 const TH1* shapeCorr,
928 Double_t trigEff,
929 Bool_t symmetrice,
930 Int_t rebin,
931 Bool_t corrEmpty,
932 Bool_t cutEdges,
933 Double_t vzMin,
934 Double_t vzMax,
935 Int_t triggerMask)
936{
937 //
938 // End of processing
939 //
940 // Parameters:
941 // sums List of sums
942 // results Output list of results
943 // shapeCorr Shape correction or nil
944 // trigEff Trigger efficiency
945 // symmetrice Whether to symmetrice the results
946 // rebin Whether to rebin the results
947 // corrEmpty Whether to correct for empty bins
948 // cutEdges Whether to cut edges when rebinning
949 // vzMin Minimum IP z coordinate
950 // vzMax Maximum IP z coordinate
951 // triggerMask Trigger mask
952 //
953
954 fSums = dynamic_cast<TList*>(sums->FindObject(GetListName()));
955 if(!fSums) {
956 AliError("Could not retrieve TList fSums");
957 return;
958 }
959
960 fOutput = new TList;
961 fOutput->SetName(GetListName());
962 fOutput->SetOwner();
963 results->Add(fOutput);
964
965 fSum = static_cast<TH2D*>(fSums->FindObject(GetName()));
966 fSumMC = static_cast<TH2D*>(fSums->FindObject(Form("%sMC", GetName())));
967 fTriggers = static_cast<TH1D*>(fSums->FindObject("triggers"));
968
969 if (!fTriggers) {
970 AliError("Couldn't find histogram 'triggers' in list");
971 return;
972 }
973 if (!fSum) {
e308a636 974 AliError(Form("Couldn't find histogram '%s' in list", GetName()));
e1f47419 975 return;
976 }
977
978 Int_t nAll = Int_t(fTriggers->GetBinContent(kAll));
979 Int_t nB = Int_t(fTriggers->GetBinContent(kB));
980 Int_t nA = Int_t(fTriggers->GetBinContent(kA));
981 Int_t nC = Int_t(fTriggers->GetBinContent(kC));
982 Int_t nE = Int_t(fTriggers->GetBinContent(kE));
983 Int_t nMB = Int_t(fTriggers->GetBinContent(kMB));
984 Int_t nTriggered = Int_t(fTriggers->GetBinContent(kWithTrigger));
985 Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
986 Int_t nAccepted = Int_t(fTriggers->GetBinContent(kAccepted));
987 Int_t nGood = nB - nA - nC + 2 * nE;
e308a636 988
e1f47419 989
e308a636 990 if (nTriggered <= 0.1) {
e1f47419 991 AliError("Number of triggered events <= 0");
992 return;
993 }
e308a636 994 if (nWithVertex <= 0.1) {
995 AliError("Number of events with vertex <= 0");
996 return;
997 }
998 if (nGood <= 0.1) {
e1f47419 999 AliWarning(Form("Number of good events=%d=%d-%d-%d+2*%d<=0",
1000 nGood, nB, nA, nC, nE));
1001 nGood = nMB;
1002 }
e308a636 1003
1004
1005 // Scaling
1006 // N_A + N_A/N_V (N_T - N_V) = N_A + N_A/N_V*N_T - N_A
1007 // = N_A/N_V * N_T
1008 //
1009 // where
1010 // N_A = nAccepted
1011 // N_V = nWithVertex
1012 // N_T = nTriggered
1013 // so that the normalisation is simply
1014 // N_A / E_V
1015 // where
1016 // E_V=N_V/N_T is the vertex efficiency from data
1017 //
1018 // Double_t alpha = Double_t(nAccepted) / nWithVertex;
1019 // Double_t vNorm = nAccepted + alpha*(nTriggered - nWithVertex);
1020 AliInfo(Form("Calculating event normalisation as "
1021 "1 / trigEff * nAccepted * nTriggered / nWithVertex "
1022 "= 1/%f * %d * %d / %d = %f",
1023 trigEff, nAccepted, nTriggered, nWithVertex,
1024 1. / trigEff * nAccepted * double(nTriggered) / nWithVertex));
1025 Double_t vNorm = double(nWithVertex) / double(nTriggered);
1026 Double_t vtxEff = vNorm;
1027 Double_t ntotal = vtxEff * trigEff;
1028
1029
e1f47419 1030 AliInfo(Form("Total of %9d events for %s\n"
e308a636 1031 " of these %9d are triggered minimum bias\n"
e1f47419 1032 " of these %9d has a %s trigger\n"
1033 " of these %9d has a vertex\n"
1034 " of these %9d were in [%+4.1f,%+4.1f]cm\n"
1035 " Triggers by type:\n"
1036 " B = %9d\n"
1037 " A|C = %9d (%9d+%-9d)\n"
1038 " E = %9d\n"
1039 " Implies %9d good triggers\n"
e308a636 1040 " Vertex efficiency: %f\n"
1041 " Trigger efficiency: %f\n"
1042 " Total number of events: %f\n",
e1f47419 1043 nAll, GetTitle(), nMB, nTriggered,
1044 AliAODForwardMult::GetTriggerString(triggerMask),
1045 nWithVertex, nAccepted,
1046 vzMin, vzMax,
e308a636 1047 nB, nA+nC, nA, nC, nE, nGood, vtxEff, trigEff, ntotal));
e1f47419 1048
1049 const char* name = GetName();
1050
1051 // Get acceptance normalisation from underflow bins
1052 TH1D* norm = ProjectX(fSum, Form("norm%s",name), 0, 0, corrEmpty, false);
1053 norm->SetDirectory(0);
1054
1055 // Scale by shape correction
1056 if(shapeCorr) fSum->Divide(shapeCorr);
1057 else AliInfo("No shape correction specified, or disabled");
1058
1059 // Project on X axis
1060 TH1D* dndeta = ProjectX(fSum, Form("dndeta%s",name), 1, fSum->GetNbinsY(),
1061 corrEmpty);
1062 dndeta->SetDirectory(0);
1063
1064 // Normalize to the acceptance -
1065 dndeta->Divide(norm);
1066
1067 // Scale by the vertex efficiency
e308a636 1068 dndeta->Scale(ntotal, "width");
e1f47419 1069
1070 SetHistogramAttributes(dndeta, kRed+1, 20, Form("ALICE %s", name));
1071 SetHistogramAttributes(norm, kRed+1,20,Form("ALICE %s normalisation", name));
1072
1073 fOutput->Add(fTriggers->Clone());
1074 if (symmetrice) fOutput->Add(Symmetrice(dndeta));
1075 fOutput->Add(dndeta);
1076 fOutput->Add(norm);
1077 fOutput->Add(Rebin(dndeta, rebin, cutEdges));
1078 if (symmetrice) fOutput->Add(Symmetrice(Rebin(dndeta, rebin, cutEdges)));
1079
1080 if (fSumMC) {
1081 // Get acceptance normalisation from underflow bins
1082 TH1D* normMC = ProjectX(fSumMC,Form("norm%sMC", name), 0, 1,
1083 corrEmpty, false);
1084 // Project onto eta axis - _ignoring_underflow_bins_!
1085 TH1D* dndetaMC = ProjectX(fSumMC,Form("dndeta%sMC", name),1,
1086 fSum->GetNbinsY(), corrEmpty);
1087 // Normalize to the acceptance
1088 dndetaMC->Divide(normMC);
1089 // Scale by the vertex efficiency
1090 dndetaMC->Scale(vNorm, "width");
1091 normMC->Scale(1. / nAccepted);
1092
1093 SetHistogramAttributes(dndetaMC, kRed+3, 21, Form("ALICE %s (MC)",name));
1094 SetHistogramAttributes(normMC, kRed+3, 21, Form("ALICE %s (MC) "
1095 "normalisation",name));
1096
1097 fOutput->Add(dndetaMC);
1098 if (symmetrice) fOutput->Add(Symmetrice(dndetaMC));
1099
1100 fOutput->Add(normMC);
1101 fOutput->Add(Rebin(dndetaMC, rebin, cutEdges));
1102
1103 if (symmetrice)
1104 fOutput->Add(Symmetrice(Rebin(dndetaMC, rebin, cutEdges)));
1105 }
e308a636 1106 if (!IsAllBin()) return;
1107
1108 // Temporary stuff
1109 TFile* forward = TFile::Open("forward.root", "READ");
1110 if (!forward) {
1111 AliWarning(Form("No forward.root file found"));
1112 return;
1113 }
1114
1115 TH1D* shapeCorrProj = 0;
1116 if (shapeCorr) {
1117 shapeCorrProj = static_cast<const TH2D*>(shapeCorr)->ProjectionX();
1118 shapeCorrProj->Scale(1. / shapeCorr->GetNbinsY());
1119 shapeCorrProj->SetDirectory(0);
1120 fOutput->Add(shapeCorrProj);
1121 }
1122
1123 TList* official = static_cast<TList*>(forward->Get("official"));
1124 if (official) {
1125 TH1F* histEta = static_cast<TH1F*>(official->FindObject("fHistEta"));
1126 if (histEta) {
1127 TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
1128 histEta->GetXaxis()->GetXmin(),
1129 histEta->GetXaxis()->GetXmax());
1130 for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
1131 oEta->SetBinContent(i, histEta->GetBinContent(i));
1132 oEta->SetBinError(i, histEta->GetBinError(i));
1133 }
1134 if (shapeCorrProj) oEta->Divide(shapeCorrProj);
1135 oEta->Scale(ntotal/nAccepted, "width");
1136 oEta->SetDirectory(0);
1137 oEta->SetMarkerStyle(21);
1138 oEta->SetMarkerColor(kCyan+3);
1139 fOutput->Add(oEta);
1140 fOutput->Add(Rebin(oEta, rebin, false));
1141 }
1142 else
1143 AliWarning(Form("Couldn't find histogram fHistEta in list %s",
1144 official->GetName()));
1145 }
1146 else
1147 AliWarning(Form("Couldn't find list 'official' in %s",forward->GetName()));
1148
1149 TList* tracks = static_cast<TList*>(forward->Get("tracks"));
1150 if (tracks) {
1151 TH1F* histEta = static_cast<TH1F*>(tracks->FindObject("fHistEta"));
1152 if (histEta) {
1153 TH1D* oEta = new TH1D("tracks", "", histEta->GetNbinsX(),
1154 histEta->GetXaxis()->GetXmin(),
1155 histEta->GetXaxis()->GetXmax());
1156 for (Int_t i = 1; i < histEta->GetNbinsX(); i++) {
1157 oEta->SetBinContent(i, histEta->GetBinContent(i));
1158 oEta->SetBinError(i, histEta->GetBinError(i));
1159 }
1160 if (shapeCorrProj) oEta->Divide(shapeCorrProj);
1161 oEta->Scale(ntotal/nAccepted, "width");
1162 oEta->SetDirectory(0);
1163 oEta->SetMarkerStyle(22);
1164 oEta->SetMarkerColor(kMagenta+2);
1165 fOutput->Add(oEta);
1166 fOutput->Add(Rebin(oEta, rebin, false));
1167 }
1168 else
1169 AliWarning(Form("Couldn't find histogram fHistEta in list %s",
1170 tracks->GetName()));
1171 }
1172 else
1173 AliWarning(Form("Couldn't find list 'tracks' in %s",forward->GetName()));
1174
1175 forward->Close();
e1f47419 1176}
1177
fb3430ac 1178//
1179// EOF
1180//