]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/FORWARD/analysis2/AliFMDHistCollector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added scripts for Casper/Valentina style P(N_{ch}) analysis
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDHistCollector.cxx
CommitLineData
7984e5f7 1//
2// This class collects the event histograms into single histograms,
3// one for each ring in each vertex bin.
4//
5// Input:
6// - AliESDFMD object possibly corrected for sharing
7//
8// Output:
9// - 5 RingHistos objects - each with a number of vertex dependent
10// 2D histograms of the inclusive charge particle density
11//
12// HistCollector used:
13// - AliFMDCorrSecondaryMap
14//
7e4038b5 15#include "AliFMDHistCollector.h"
16#include <AliESDFMD.h>
17#include <TAxis.h>
18#include <TList.h>
19#include <TMath.h>
0bd4b00f 20#include "AliForwardCorrectionManager.h"
8449e3e0 21#include "AliFMDCorrSecondaryMap.h"
7e4038b5 22#include "AliLog.h"
23#include <TH2D.h>
4f9319f3 24#include <TH3D.h>
7e4038b5 25#include <TH1I.h>
26#include <TProfile.h>
27#include <TProfile2D.h>
8449e3e0 28#include <TObjArray.h>
7e4038b5 29#include <TArrayI.h>
0bd4b00f 30#include <TROOT.h>
31#include <iostream>
32#include <iomanip>
7e4038b5 33
34ClassImp(AliFMDHistCollector)
35#if 0
36; // For Emacs
37#endif
38
5bb5d1f6 39//____________________________________________________________________
40AliFMDHistCollector::AliFMDHistCollector()
41 : fNCutBins(0),
42 fCorrectionCut(0),
5bb5d1f6 43 fDebug(0),
44 fList(0),
45 fSumRings(0),
46 fCoverage(0),
47 fMergeMethod(kStraightMean),
ce731b9f 48 fFiducialMethod(kByCut),
49 fSkipFMDRings(0),
8449e3e0 50 fBgAndHitMaps(false),
51 fVtxList(0),
52 fByCent(0),
53 fDoByCent(false)
6ab100ec 54{
5ca83fee 55 DGUARD(fDebug, 3, "Default CTOR of AliFMDHistCollector");
6ab100ec 56}
5bb5d1f6 57
58//____________________________________________________________________
59AliFMDHistCollector::AliFMDHistCollector(const char* title)
60 : TNamed("fmdHistCollector", title),
61 fNCutBins(2),
62 fCorrectionCut(0.5),
5bb5d1f6 63 fDebug(0),
64 fList(0),
65 fSumRings(0),
66 fCoverage(0),
67 fMergeMethod(kStraightMean),
ce731b9f 68 fFiducialMethod(kByCut),
69 fSkipFMDRings(0),
8449e3e0 70 fBgAndHitMaps(false),
71 fVtxList(0),
72 fByCent(0),
73 fDoByCent(false)
5bb5d1f6 74{
5ca83fee 75 DGUARD(fDebug, 3, "Named CTOR of AliFMDHistCollector: %s", title);
5bb5d1f6 76}
77//____________________________________________________________________
78AliFMDHistCollector::AliFMDHistCollector(const AliFMDHistCollector& o)
79 : TNamed(o),
80 fNCutBins(o.fNCutBins),
81 fCorrectionCut(o.fCorrectionCut),
5bb5d1f6 82 fDebug(o.fDebug),
83 fList(o.fList),
84 fSumRings(o.fSumRings),
85 fCoverage(o.fCoverage),
86 fMergeMethod(o.fMergeMethod),
ce731b9f 87 fFiducialMethod(o.fFiducialMethod),
88 fSkipFMDRings(o.fSkipFMDRings),
8449e3e0 89 fBgAndHitMaps(o.fBgAndHitMaps),
90 fVtxList(o.fVtxList),
91 fByCent(o.fByCent),
92 fDoByCent(o.fDoByCent)
6ab100ec 93{
5ca83fee 94 DGUARD(fDebug, 3, "Copy CTOR of AliFMDHistCollector");
6ab100ec 95}
5bb5d1f6 96
12fffad7 97//____________________________________________________________________
98AliFMDHistCollector::~AliFMDHistCollector()
99{
6ab100ec 100 DGUARD(fDebug, 3, "DTOR of AliFMDHistCollector");
b7ab8a2c 101 // if (fList) delete fList;
12fffad7 102}
7e4038b5 103//____________________________________________________________________
104AliFMDHistCollector&
105AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
106{
7984e5f7 107 //
108 // Assignement operator
109 //
110 // Parameters:
111 // o Object to assign from
112 //
113 // Return:
114 // Reference to this object
115 //
6ab100ec 116 DGUARD(fDebug, 3, "Assignment of AliFMDHistCollector");
d015ecfe 117 if (&o == this) return *this;
ea3e5d95 118 TNamed::operator=(o);
7e4038b5 119
120 fNCutBins = o.fNCutBins;
121 fCorrectionCut = o.fCorrectionCut;
ea3e5d95 122 fDebug = o.fDebug;
12fffad7 123 fList = o.fList;
124 fSumRings = o.fSumRings;
125 fCoverage = o.fCoverage;
e1f47419 126 fMergeMethod = o.fMergeMethod;
127 fFiducialMethod = o.fFiducialMethod;
ce731b9f 128 fSkipFMDRings = o.fSkipFMDRings;
129 fBgAndHitMaps = o.fBgAndHitMaps;
8449e3e0 130 fVtxList = o.fVtxList;
131 fByCent = o.fByCent;
132 fDoByCent = o.fDoByCent;
7e4038b5 133 return *this;
134}
135
136//____________________________________________________________________
137void
5934a3e3 138AliFMDHistCollector::SetupForData(const TAxis& vtxAxis,
8449e3e0 139 const TAxis& etaAxis)
7e4038b5 140{
7984e5f7 141 //
142 // Intialise
143 //
144 // Parameters:
145 // vtxAxis Vertex axis
146 //
6ab100ec 147 DGUARD(fDebug, 1, "Initialization of AliFMDHistCollector");
7984e5f7 148
8449e3e0 149 // AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5 150
12fffad7 151 fSumRings = new TH2D("sumRings", "Sum in individual rings",
152 etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
153 5, 1, 6);
154 fSumRings->Sumw2();
155 fSumRings->SetDirectory(0);
156 fSumRings->SetXTitle("#eta");
157 fSumRings->GetYaxis()->SetBinLabel(1,"FMD1i");
158 fSumRings->GetYaxis()->SetBinLabel(2,"FMD2i");
159 fSumRings->GetYaxis()->SetBinLabel(3,"FMD2o");
160 fSumRings->GetYaxis()->SetBinLabel(4,"FMD3i");
161 fSumRings->GetYaxis()->SetBinLabel(5,"FMD3o");
162 fList->Add(fSumRings);
163
164 fCoverage = new TH2D("coverage", "#eta coverage per v_{z}",
e1f47419 165 etaAxis.GetNbins(),etaAxis.GetXmin(),etaAxis.GetXmax(),
166 vtxAxis.GetNbins(),vtxAxis.GetXmin(),vtxAxis.GetXmax());
12fffad7 167 fCoverage->SetDirectory(0);
168 fCoverage->SetXTitle("#eta");
169 fCoverage->SetYTitle("v_{z} [cm]");
170 fCoverage->SetZTitle("n_{bins}");
171 fList->Add(fCoverage);
172
7e4038b5 173
6af747ad 174 // --- Add parameters to output ------------------------------------
175 fList->Add(AliForwardUtil::MakeParameter("nCutBins",fNCutBins));
176 fList->Add(AliForwardUtil::MakeParameter("skipRings",fSkipFMDRings));
177 fList->Add(AliForwardUtil::MakeParameter("bgAndHits",fBgAndHitMaps));
178 fList->Add(AliForwardUtil::MakeParameter("fiducial",Int_t(fFiducialMethod)));
179 fList->Add(AliForwardUtil::MakeParameter("fiducialCut",fCorrectionCut));
180 fList->Add(AliForwardUtil::MakeParameter("merge",Int_t(fMergeMethod)));
181
8449e3e0 182 UShort_t nVz = vtxAxis.GetNbins();
183 fVtxList = new TObjArray(nVz, 1);
184 fVtxList->SetName("histCollectorVtxBins");
185 fVtxList->SetOwner();
186
7e4038b5 187 // Find the eta bin ranges
0bd4b00f 188 for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
e1f47419 189 Double_t vMin = vtxAxis.GetBinLowEdge(iVz);
190 Double_t vMax = vtxAxis.GetBinUpEdge(iVz);
8449e3e0 191 VtxBin* bin = new VtxBin(iVz, vMin, vMax, fNCutBins);
192 fVtxList->AddAt(bin, iVz);
7e4038b5 193
8449e3e0 194 bin->SetupForData(fCoverage, fSkipFMDRings, fFiducialMethod,
195 fCorrectionCut, fList, etaAxis,
196 fBgAndHitMaps, fBgAndHitMaps);
197 }
7e4038b5 198
8449e3e0 199 if (!fDoByCent) return;
200
201 fByCent = new TList;
202 fByCent->SetName("byCentrality");
203 fByCent->SetOwner();
204 fList->Add(fByCent);
205
206 Int_t nCent = 101;
207 Double_t minCent = -.5;
208 Double_t maxCent = 100.5;
209 for (Int_t i = 0; i < 5; i++) {
210 UShort_t d;
211 Char_t r;
212 GetDetRing(i, d, r);
213
214 TH3* h = new TH3D(Form("FMD%d%c", d, r),
215 Form("dN/d#eta per centrality for FMD%d%c", d, r),
216 etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(),
217 nCent, minCent, maxCent, 1, 0, 1);
218 h->SetXTitle("#eta");
219 h->SetYTitle("Centrality [%]");
220 h->SetZTitle("dN/d#eta");
221 h->SetDirectory(0);
222 h->SetMarkerColor(AliForwardUtil::RingColor(d, r));
223 h->SetMarkerStyle(20);
224 fByCent->Add(h);
7e4038b5 225 }
226}
e1f47419 227//____________________________________________________________________
228Bool_t
8449e3e0 229AliFMDHistCollector::CheckCorrection(FiducialMethod m,
230 Double_t cut,
231 const TH2D* bg,
232 Int_t ie,
233 Int_t ip)
e1f47419 234{
235 //
236 // Check if we should include the bin in the data range
237 //
238 // Parameters:
239 // bg Secondary map histogram
240 // ie Eta bin
241 // ip Phi bin
242 //
243 // Return:
244 // True if to be used
245 //
246 Double_t c = bg->GetBinContent(ie,ip);
8449e3e0 247 switch (m) {
e1f47419 248 case kByCut:
8449e3e0 249 return c >= cut;
e1f47419 250 case kDistance:
251 if (2 * c < bg->GetBinContent(ie+1,ip) ||
252 2 * c < bg->GetBinContent(ie-1,ip)) return false;
253 return true;
254 default:
8449e3e0 255 AliErrorClass("No fiducal cut method defined");
e1f47419 256 }
257 return false;
258}
259
12fffad7 260//____________________________________________________________________
261void
5934a3e3 262AliFMDHistCollector::CreateOutputObjects(TList* dir)
12fffad7 263{
264 //
265 // Output diagnostic histograms to directory
266 //
267 // Parameters:
268 // dir List to write in
269 //
6ab100ec 270 DGUARD(fDebug, 1, "Define output of AliFMDHistCollector");
12fffad7 271 fList = new TList;
272 fList->SetOwner();
273 fList->SetName(GetName());
274 dir->Add(fList);
5bb5d1f6 275
12fffad7 276}
277
8449e3e0 278//____________________________________________________________________
279Bool_t
280AliFMDHistCollector::CheckSkip(UShort_t d, Char_t r, UShort_t skips)
281{
8449e3e0 282 UShort_t q = (r == 'I' || r == 'i' ? 0 : 1);
bfab35d9 283 UShort_t c = 1 << (d-1);
284 UShort_t t = 1 << (c+q-1);
285
8449e3e0 286 return (t & skips) == t;
287}
5bb5d1f6 288
7e4038b5 289//____________________________________________________________________
290Int_t
8449e3e0 291AliFMDHistCollector::GetIdx(UShort_t d, Char_t r)
7e4038b5 292{
7984e5f7 293 //
294 // Get the ring index from detector number and ring identifier
295 //
296 // Parameters:
297 // d Detector
298 // r Ring identifier
299 //
300 // Return:
301 // ring index or -1 in case of problems
302 //
7e4038b5 303 Int_t idx = -1;
304 switch (d) {
305 case 1: idx = 0; break;
306 case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
307 case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
308 }
309 return idx;
310}
311//____________________________________________________________________
312void
8449e3e0 313AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r)
7e4038b5 314{
7984e5f7 315 //
316 // Get the detector and ring from the ring index
317 //
318 // Parameters:
319 // idx Ring index
320 // d On return, the detector or 0 in case of errors
321 // r On return, the ring id or '0' in case of errors
322 //
7e4038b5 323 d = 0;
324 r = '\0';
325 switch (idx) {
326 case 0: d = 1; r = 'I'; break;
327 case 1: d = 2; r = 'I'; break;
328 case 2: d = 2; r = 'O'; break;
329 case 3: d = 3; r = 'I'; break;
330 case 4: d = 3; r = 'O'; break;
331 }
332}
333
334//____________________________________________________________________
8449e3e0 335AliFMDHistCollector::VtxBin*
336AliFMDHistCollector::GetVtxBin(Int_t ivtx)
7e4038b5 337{
7984e5f7 338 // Parameters:
7984e5f7 339 // vtxBin Vertex bin (1 based)
8449e3e0 340 if (!fVtxList) return 0;
341 if (ivtx < 1 || ivtx > fVtxList->GetEntriesFast()) return 0;
342 VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(ivtx));
343 return bin;
7e4038b5 344}
7e4038b5 345//____________________________________________________________________
8449e3e0 346const AliFMDHistCollector::VtxBin*
347AliFMDHistCollector::GetVtxBin(Int_t ivtx) const
7e4038b5 348{
7984e5f7 349 // Parameters:
8449e3e0 350 // vtxBin Vertex bin (1 based)
351 if (!fVtxList) return 0;
352 if (ivtx < 1 || ivtx > fVtxList->GetEntriesFast()) return 0;
353 VtxBin* bin = static_cast<VtxBin*>(fVtxList->At(ivtx));
354 return bin;
7e4038b5 355}
7e4038b5 356
e1f47419 357//____________________________________________________________________
358void
8449e3e0 359AliFMDHistCollector::MergeBins(MergeMethod m,
360 Double_t c, Double_t e,
e1f47419 361 Double_t oc, Double_t oe,
8449e3e0 362 Double_t& rc, Double_t& re)
e1f47419 363{
364 //
365 // Merge bins accoring to set method
366 //
367 // Parameters:
368 // c Current content
369 // e Current error
370 // oc Old content
371 // oe Old error
372 // rc On return, the new content
373 // re On return, tne new error
374 //
375 rc = re = 0;
8449e3e0 376 switch (m) {
e1f47419 377 case kStraightMean:
bfab35d9 378 case kPreferInner: // We only get these two when there's an overlap
379 case kPreferOuter: // between like rings, and we should do a straight mean
e1f47419 380 // calculate the average of old value (half the original),
381 // and this value, as well as the summed squared errors
382 // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
383 // and half the error of this.
384 //
385 // So, on the first overlapping histogram we get
386 //
387 // c = c_1 / 2
388 // e = sqrt((e_1 / 2)^2) = e_1/2
389 //
390 // On the second we get
391 //
392 // c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
393 // e' = sqrt(e^2 + (e_2/2)^2)
394 // = sqrt(e_1^2/4 + e_2^2/4)
395 // = sqrt(1/4 * (e_1^2+e_2^2))
396 // = 1/2 * sqrt(e_1^2 + e_2^2)
397 rc = oc + c/2;
398 re = TMath::Sqrt(oe*oe+(e*e)/4);
399 break;
400 case kStraightMeanNoZero:
401 // If there's data in the overlapping histogram,
402 // calculate the average and add the errors in
403 // quadrature.
404 //
405 // If there's no data in the overlapping histogram,
406 // then just fill in the data
407 if (oe > 0) {
408 rc = (oc + c)/2;
409 re = TMath::Sqrt(oe*oe + e*e)/2;
410 }
411 else {
412 rc = c;
413 re = e;
414 }
415 break;
416 case kWeightedMean: {
417 // Calculate the weighted mean
418 Double_t w = 1/(e*e);
419 Double_t sc = w * c;
420 Double_t sw = w;
421 if (oe > 0) {
422 Double_t ow = 1/(oe*oe);
423 sc += ow * oc;
424 sw += ow;
425 }
426 rc = sc / sw;
427 re = TMath::Sqrt(1 / sw);
428 }
429 break;
430 case kLeastError:
431 if (e < oe) {
432 rc = c;
433 re = e;
434 }
435 else {
436 rc = oc;
437 re = oe;
438 }
439 break;
5ca83fee 440 case kSum:
441 rc = c + oc;
442 re = TMath::Sqrt(oe * oe + e * e);//Add in quadarature
443 break;
e1f47419 444 default:
8449e3e0 445 AliErrorClass("No method for defining content of overlapping bins defined");
e1f47419 446 return;
447 }
448}
449
7e4038b5 450//____________________________________________________________________
451Bool_t
5bb5d1f6 452AliFMDHistCollector::Collect(const AliForwardUtil::Histos& hists,
453 AliForwardUtil::Histos& sums,
0bd4b00f 454 UShort_t vtxbin,
4f9319f3 455 TH2D& out,
8449e3e0 456 Double_t cent)
7e4038b5 457{
7984e5f7 458 //
459 // Do the calculations
460 //
461 // Parameters:
462 // hists Cache of histograms
463 // vtxBin Vertex bin (1 based)
464 // out Output histogram
465 //
466 // Return:
467 // true on successs
468 //
6ab100ec 469 DGUARD(fDebug, 1, "Collect final histogram of AliFMDHistCollector");
8449e3e0 470 // AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
471 // const TAxis* vtxAxis = fcm.GetVertexAxis();
472 // Double_t vMin = vtxAxis->GetBinLowEdge(vtxbin);
473 // Double_t vMax = vtxAxis->GetBinUpEdge(vtxbin);
474 VtxBin* bin = GetVtxBin(vtxbin);
717a4c2c 475 if (!bin) return false;
8449e3e0 476 Bool_t ret = bin->Collect(hists, sums, out, fSumRings, cent,
477 fMergeMethod, fSkipFMDRings,
478 fByCent);
479
480 return ret;
481}
482
483//____________________________________________________________________
484void
485AliFMDHistCollector::Print(Option_t* /* option */) const
486{
487 //
488 // Print information
489 //
490 // Parameters:
491 // option Not used
492 //
493 char ind[gROOT->GetDirLevel()+1];
494 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
495 ind[gROOT->GetDirLevel()] = '\0';
496 std::cout << ind << ClassName() << ": " << GetName() << '\n'
497 << ind << " # of cut bins: " << fNCutBins << '\n'
498 << ind << " Fiducal method: "
499 << (fFiducialMethod == kByCut ? "cut" : "distance") << "\n"
500 << ind << " Fiducial cut: " << fCorrectionCut << "\n"
501 << ind << " Merge method: ";
502 switch (fMergeMethod) {
503 case kStraightMean: std::cout << "straight mean\n"; break;
504 case kStraightMeanNoZero: std::cout << "straight mean (no zeros)\n"; break;
505 case kWeightedMean: std::cout << "weighted mean\n"; break;
506 case kLeastError: std::cout << "least error\n"; break;
507 case kSum: std::cout << "straight sum\n"; break;
bfab35d9 508 case kPreferInner: std::cout << "prefer inners\n"; break;
509 case kPreferOuter: std::cout << "prefer outers\n"; break;
8449e3e0 510 }
511
512 if (!fVtxList) return;
513
514 std::cout << ind << " Bin ranges:\n" << ind << " rings | Range ";
515 Int_t nVz = fVtxList->GetEntriesFast();
516 for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
517 UShort_t d = 0;
518 Char_t r = 0;
519 GetDetRing(iIdx, d, r);
520 std::cout << ind << " | FMD" << d << r << " ";
521 }
522 std::cout << '\n' << ind << " /vz_bin |-----------";
523 for (Int_t iIdx = 0; iIdx < 5; iIdx++)
524 std::cout << "-+--------";
525 std::cout << std::endl;
526
527 for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
528 const VtxBin* bin = GetVtxBin(iVz);
529 if (!bin) continue;
530 std::cout << " " << std::right << std::setw(6) << iVz << " | "
531 << std::setw(3) << bin->fLow << " - " << std::left
532 << std::setw(3) << bin->fHigh << " ";
533 for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
534 Int_t first, last;
535 bin->GetFirstAndLast(iIdx, first, last);
536 std::cout << " | " << std::setw(3) << first << "-"
537 << std::setw(3) << last;
538 }
539 std::cout << std::endl;
540 }
541}
542
543//____________________________________________________________________
544AliFMDHistCollector::VtxBin::VtxBin(Int_t idx, Double_t minIpZ, Double_t maxIpZ,
545 Int_t nCutBins)
546 : fIndex(idx),
547 fLow(minIpZ),
548 fHigh(maxIpZ),
549 fHitMap(0),
550 fFirstBin(1),
551 fLastBin(1),
552 fNCutBins(nCutBins)
553{
554}
555//____________________________________________________________________
556AliFMDHistCollector::VtxBin::VtxBin(const VtxBin& o)
557 : TObject(o),
558 fIndex(o.fIndex),
559 fLow(o.fLow),
560 fHigh(o.fHigh),
561 fHitMap(o.fHitMap),
562 fFirstBin(o.fFirstBin),
563 fLastBin(o.fLastBin),
564 fNCutBins(o.fNCutBins)
565{
566}
567//____________________________________________________________________
568AliFMDHistCollector::VtxBin&
569AliFMDHistCollector::VtxBin::operator=(const VtxBin& o)
570{
571 if (&o == this) return *this;
572 fIndex = o.fIndex;
573 fLow = o.fLow;
574 fHigh = o.fHigh;
575 fHitMap = o.fHitMap;
576 fFirstBin = o.fFirstBin;
577 fLastBin = o.fLastBin;
578 fNCutBins = o.fNCutBins;
579 return *this;
580}
581//____________________________________________________________________
582const Char_t*
583AliFMDHistCollector::VtxBin::GetName() const
584{
585 return Form("%c%03d_%c%03d",
586 (fLow >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fLow)),
587 (fHigh >= 0 ? 'p' : 'm'), Int_t(TMath::Abs(fHigh)));
588}
589//____________________________________________________________________
590void
591AliFMDHistCollector::VtxBin::SetupForData(TH2* coverage,
592 UShort_t skips,
593 FiducialMethod fiducial,
594 Double_t cut,
595 TList* l,
596 const TAxis& etaAxis,
597 Bool_t doHitMaps,
598 Bool_t storeSecMap)
599{
600 TList* out = 0;
601 if (doHitMaps || storeSecMap) {
602 out = new TList;
603 out->SetName(GetName());
604 out->SetOwner();
605 l->Add(out);
606 }
607 if (doHitMaps) {
608 fHitMap = new AliForwardUtil::Histos();
609 fHitMap->Init(etaAxis);
610 }
611 fFirstBin.Set(5);
612 fLastBin.Set(5);
613
73e536c6 614 AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
73e536c6 615
8449e3e0 616 for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
617 UShort_t d = 0;
618 Char_t r = 0;
619 GetDetRing(iIdx, d, r);
620
621 // Skipping selected FMD rings
622 if (CheckSkip(d, r, skips)) continue;
623
624 // Get the background object
625 TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,UShort_t(fIndex));
626 Int_t nEta = bg->GetNbinsX();
627 Int_t first = nEta+1;
628 Int_t last = 0;
629
630 // Loop over the eta bins
631 for (Int_t ie = 1; ie <= nEta; ie++) {
632 // Loop over the phi bins to make sure that we
633 // do not have holes in the coverage
634 bool ok = true;
635 for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
636 if (!CheckCorrection(fiducial, cut, bg, ie, ip)) {
637 ok = false;
638 continue;
639 }
640 }
641 if (!ok) continue;
642
643 first = TMath::Min(ie, first);
644 last = TMath::Max(ie, last);
645 }
646 // Store result of first/last bin for this ring
647 fFirstBin[iIdx] = first;
648 fLastBin[iIdx] = last;
649
650 if (fHitMap) {
651 TH2* h = fHitMap->Get(d, r);
652 h->SetDirectory(0);
653 h->SetName(Form("hitMapFMD%d%c", d, r));
654 // if (r == 'O') h->RebinY(2);
655 out->Add(h);
656 }
657
658 TH2D* obg=0;
659 if(storeSecMap) {
660 obg = static_cast<TH2D*>(bg->Clone(Form("secMapFMD%d%c", d, r)));
661 obg->SetDirectory(0);
662 obg->Reset();
663 out->Add(obg);
664 }
665 // Fill diagnostics histograms
666 for (Int_t ie = first+fNCutBins; ie <= last-fNCutBins; ie++) {
667 Double_t old = coverage->GetBinContent(ie, fIndex);
668 coverage->SetBinContent(ie, fIndex, old+1);
669 if(obg) {
670 for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
671 obg->SetBinContent(ie, ip, bg->GetBinContent(ie, ip));
672 obg->SetBinError(ie, ip, bg->GetBinError(ie, ip));
673 } // for (ip)
674 } // if (doSecHits)
675 } // for (ie)
676 } // for (iIdx)
677}
678
679//____________________________________________________________________
680void
681AliFMDHistCollector::VtxBin::GetFirstAndLast(Int_t idx,
682 Int_t& first,
683 Int_t& last) const
684{
685 // Get the first and last eta bin to use for a given ring and vertex
686 //
687 // Parameters:
688 // idx Ring index as given by GetIdx
689 // first On return, the first eta bin to use
690 // last On return, the last eta bin to use
691 //
692 first = 0;
693 last = 0;
694
695 if (idx < 0 || idx >= fFirstBin.GetSize()) return;
696
697 first = fFirstBin.At(idx)+fNCutBins;
698 last = fLastBin.At(idx)-fNCutBins;
699}
700//____________________________________________________________________
701Int_t
702AliFMDHistCollector::VtxBin::GetFirst(Int_t idx) const
703{
704 Int_t first, last;
705 GetFirstAndLast(idx, first , last);
706 return first;
707}
708//____________________________________________________________________
709Int_t
710AliFMDHistCollector::VtxBin::GetLast(Int_t idx) const
711{
712 Int_t first, last;
713 GetFirstAndLast(idx, first , last);
714 return last;
715}
716
717//____________________________________________________________________
718Bool_t
719AliFMDHistCollector::VtxBin::Collect(const AliForwardUtil::Histos& hists,
720 AliForwardUtil::Histos& sums,
721 TH2D& out,
722 TH2D* sumRings,
723 Double_t cent,
724 MergeMethod m,
725 UShort_t skips,
726 TList* byCent)
727{
7e4038b5 728 for (UShort_t d=1; d<=3; d++) {
729 UShort_t nr = (d == 1 ? 1 : 2);
730 for (UShort_t q=0; q<nr; q++) {
731 Char_t r = (q == 0 ? 'I' : 'O');
8449e3e0 732 if (CheckSkip(d, r, skips)) continue;
d9af24a9 733
7e4038b5 734 TH2D* h = hists.Get(d,r);
8449e3e0 735 TH2D* o = sums.Get(d, r);
7e4038b5 736 TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
12fffad7 737 Int_t i = (d == 1 ? 1 : 2*d + (q == 0 ? -2 : -1));
8449e3e0 738
5bb5d1f6 739 // Get valid range
7e4038b5 740 Int_t first = 0;
741 Int_t last = 0;
8449e3e0 742 GetFirstAndLast(d, r, first, last);
5bb5d1f6 743
744 // Zero outside valid range
5ca83fee 745 Int_t nY = t->GetNbinsY();
8449e3e0 746 Int_t nX = t->GetNbinsX();
5ca83fee 747 for (Int_t iPhi = 0; iPhi <= nY+1; iPhi++) {
5bb5d1f6 748 // Lower range
749 for (Int_t iEta = 1; iEta < first; iEta++) {
750 t->SetBinContent(iEta,iPhi,0);
751 t->SetBinError(iEta,iPhi,0);
752 }
8449e3e0 753 for (Int_t iEta = last+1; iEta <= nX; iEta++) {
5bb5d1f6 754 t->SetBinContent(iEta,iPhi,0);
755 t->SetBinError(iEta,iPhi,0);
756 }
757 }
bfab35d9 758 // Fill under-flow bins with eta coverage
759 for (Int_t iEta = first; iEta <= last; iEta++)
5bb5d1f6 760 t->SetBinContent(iEta,0,1);
bfab35d9 761
5bb5d1f6 762 // Add to our per-ring sum
763 o->Add(t);
bfab35d9 764
8449e3e0 765 // If we store hit maps, update here
766 if (fHitMap) fHitMap->Get(d, r)->Add(t);
4f9319f3 767
8449e3e0 768 if (byCent) {
769 TH3* dNdetaCent = static_cast<TH3*>(byCent->At(i-1));
770 if (cent >= 0 && dNdetaCent) {
771 Int_t iCent = dNdetaCent->GetYaxis()->FindBin(cent);
772
773 if (iCent > 0 && iCent <= dNdetaCent->GetNbinsY()) {
774 // Make a projection of data
775 TH1* proj = static_cast<TH1*>(t->ProjectionX("tmp", 1, nY));
776 proj->SetDirectory(0);
777 for (Int_t iEta = 1; iEta <= nX; iEta++) {
778 Double_t v1 = proj->GetBinContent(iEta);
779 Double_t e1 = proj->GetBinError(iEta);
780 Double_t v2 = dNdetaCent->GetBinContent(iEta, iCent, 1);
781 Double_t e2 = dNdetaCent->GetBinError(iEta, iCent, 1);
782 dNdetaCent->SetBinContent(iEta,iCent,1, v1+v2);
783 dNdetaCent->SetBinError(iEta,iCent,1, TMath::Sqrt(e1*e1+e2*e2));
784
785 // Check under-/overflow bins
786 Double_t uF = t->GetBinContent(iEta, 0);
787 Double_t oF = t->GetBinContent(iEta, nY+1);
788 if (uF > 0) {
789 Double_t old = dNdetaCent->GetBinContent(iEta, iCent, 0);
790 dNdetaCent->SetBinContent(iEta, iCent, 0, old + uF);
791 }
792 if (oF > 0) {
793 Double_t old = dNdetaCent->GetBinContent(iEta, iCent, 2);
794 dNdetaCent->SetBinContent(iEta, iCent, 2, old + oF);
795 }
796 } // for(iEta)
797 delete proj;
798 } // if(iCent)
799 } // if (cent)
800 } // if (byCent)
4f9319f3 801
5bb5d1f6 802 // Outer rings have better phi segmentation - rebin to same as inner.
803 if (q == 1) t->RebinY(2);
7e4038b5 804
5ca83fee 805 nY = t->GetNbinsY();
bfab35d9 806
7e4038b5 807 // Now update profile output
808 for (Int_t iEta = first; iEta <= last; iEta++) {
809
810 // Get the possibly overlapping histogram
8449e3e0 811 Int_t overlap = GetOverlap(d,r,iEta);
7e4038b5 812
5ca83fee 813 // Get factor
bfab35d9 814 MergeMethod mm = m; // Possibly override method locally
815 Float_t fac = 1;
816 if (m != kSum && overlap >= 0) {
817 // Default is to average
818 fac = 0.5;
819 if (m == kPreferInner) {
820 // Current one is an outer overlapping an inner
821 if ((r == 'o' || r == 'O') &&
822 (overlap == 0 || overlap == 1 || overlap == 3))
823 // Do not use this signal
824 fac = 0;
825 // Current one is inner overlapping an outer
826 else if ((r == 'i' || r == 'I') && (overlap == 2 || overlap == 4))
827 // Prefer this one
828 fac = 1;
829 else
830 // In case of two overlapping inners
831 mm = kStraightMean;
832 }
833 else if (m == kPreferOuter) {
834 // Current one is an inner overlapping an outer
835 if ((r == 'i' || r == 'I') && (overlap == 2 || overlap == 4))
836 // Do not use this signal
837 fac = 0;
838 else if ((r == 'O' || r == 'o') &&
839 (overlap == 0 || overlap == 1 || overlap == 3))
840 fac = 1;
841 else
842 // In case of two overlapping outers
843 mm = kStraightMean;
844 }
845 }
5ca83fee 846
847 // Fill eta acceptance for this event into the phi underflow bin
7e4038b5 848 Float_t ooc = out.GetBinContent(iEta,0);
5ca83fee 849 out.SetBinContent(iEta, 0, ooc + fac);
850
851 // Fill phi acceptance for this event into the phi overflow bin
852 Float_t oop = out.GetBinContent(iEta,nY+1);
853 Float_t nop = t->GetBinContent(iEta,nY+1);
854#if 0
855 Info("", "etaBin=%3d Setting phi acceptance to %f(%f+%f)=%f",
856 iEta, fac, oop, nop, fac*(oop+nop));
857#endif
858 out.SetBinContent(iEta, nY+1, fac * nop + oop);
7e4038b5 859
5bb5d1f6 860 // Should we loop over h or t Y bins - I think it's t
5ca83fee 861 for (Int_t iPhi = 1; iPhi <= nY; iPhi++) {
12fffad7 862 Double_t c = t->GetBinContent(iEta,iPhi);
863 Double_t e = t->GetBinError(iEta,iPhi);
864 Double_t ee = t->GetXaxis()->GetBinCenter(iEta);
8449e3e0 865 sumRings->Fill(ee, i, c);
7e4038b5 866
bfab35d9 867 // If there's no signal or the signal was ignored because we
868 // prefer the inners/outers, continue if (e <= 0) continue;
869 if (fac <= 0 || c <= 0 || e <= 0) continue;
7e4038b5 870
bfab35d9 871 // If there's no overlapping histogram (ring) or if we
872 // prefer inner/outer, then fill in data and continue to the
873 // next phi bin
874 if (overlap < 0 || fac >= 1) {
7e4038b5 875 out.SetBinContent(iEta,iPhi,c);
876 out.SetBinError(iEta,iPhi,e);
877 continue;
878 }
879
880 // Get the current bin content and error
881 Double_t oc = out.GetBinContent(iEta,iPhi);
882 Double_t oe = out.GetBinError(iEta,iPhi);
883
e1f47419 884 Double_t rc, re;
bfab35d9 885 MergeBins(mm, c, e, oc, oe, rc, re);
e1f47419 886 out.SetBinContent(iEta,iPhi, rc);
887 out.SetBinError(iEta,iPhi, re);
7e4038b5 888 }
889 }
36ffcf83 890 // Remove temporary histogram
7e4038b5 891 delete t;
892 } // for r
893 } // for d
8449e3e0 894 return true;
7e4038b5 895}
0bd4b00f 896//____________________________________________________________________
8449e3e0 897Int_t
898AliFMDHistCollector::VtxBin::GetOverlap(UShort_t d, Char_t r,
899 Int_t bin) const
0bd4b00f 900{
7984e5f7 901 //
8449e3e0 902 // Get the possibly overlapping histogram of eta bin @a e in
903 // detector and ring
7984e5f7 904 //
905 // Parameters:
8449e3e0 906 // d Detector
907 // r Ring
908 // e Eta bin
909 // v Vertex bin (1 based)
7984e5f7 910 //
8449e3e0 911 // Return:
912 // Overlapping histogram index or -1
913 //
914
915 Int_t other = -1;
916 if (d == 1) {
917 if (bin <= GetLast(2,'I')) other = GetIdx(2,'I');
e1f47419 918 }
8449e3e0 919 else if (d == 2 && r == 'I') {
920 if (bin <= GetLast(2, 'O')) other = GetIdx(2, 'O');
921 else if (bin >= GetFirst(1, 'I')) other = GetIdx(1, 'I');
5bb5d1f6 922 }
8449e3e0 923 else if (d == 2 && r == 'O') {
924 if (bin >= GetFirst(2, 'I')) other = GetIdx(2,'I');
0bd4b00f 925 }
8449e3e0 926 else if (d == 3 && r == 'O') {
927 if (bin <= GetLast(3, 'I')) other = GetIdx(3, 'I');
928 }
929 else if (d == 3 && r == 'I') {
930 if (bin >= GetFirst(3, 'O')) other = GetIdx(3, 'O');
931 }
932 // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
933 return other;
0bd4b00f 934}
8449e3e0 935
7e4038b5 936
937//____________________________________________________________________
938//
939// EOF
940//
941
942
943