]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDMCSharingFilter.cxx
CommitLineData
dc0b1641 1//
7984e5f7 2// Class to do the sharing correction for MC data.
dc0b1641 3//
7984e5f7 4// Input:
5// - AliESDFMD object - from reconstruction
6// - Kinematics
7// - Track-References
8//
9// Output:
10// - AliESDFMD object - copy of input, but with signals merged
11//
12// Corrections used:
5bb5d1f6 13// - ELoss fits
7984e5f7 14//
15// Histograms:
16// - For each ring (FMD1i, FMD2i, FMD2o, FMD3i, FMD3o) the distribution of
17// signals before and after the filter.
18// - For each ring (see above), an array of distributions of number of
19// hit strips for each vertex bin (if enabled - see Init method)
20//
dc0b1641 21#include "AliFMDMCSharingFilter.h"
22#include <AliESDFMD.h>
23#include <AliMCEvent.h>
24#include <AliTrackReference.h>
25#include <AliStack.h>
26#include <TAxis.h>
27#include <TList.h>
28#include <TH1.h>
29#include <TMath.h>
dc0b1641 30#include "AliFMDStripIndex.h"
5bb5d1f6 31#include "AliFMDFloatMap.h"
dc0b1641 32#include <AliLog.h>
33#include <TROOT.h>
34#include <iostream>
35#include <iomanip>
36
37ClassImp(AliFMDMCSharingFilter)
38#if 0
39; // This is for Emacs
40#endif
41
42
dc0b1641 43//____________________________________________________________________
44AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
45 : AliFMDSharingFilter(title),
46 fFMD1i(0),
47 fFMD2i(0),
48 fFMD2o(0),
49 fFMD3i(0),
50 fFMD3o(0),
5bb5d1f6 51 fSumEta(0),
52 fOperComp(0),
53 fThetaVsNr(0),
54 fOnlyPrimary(false)
dc0b1641 55{
7984e5f7 56 //
57 // Constructor
58 //
59 // Parameters:
60 // title Title of object - not significant
61 //
1174780f 62 fFMD1i = new TH2D("FMD1i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
63 fFMD2i = new TH2D("FMD2i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
64 fFMD2o = new TH2D("FMD2o_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
65 fFMD3i = new TH2D("FMD3i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
66 fFMD3o = new TH2D("FMD3o_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
dc0b1641 67 fFMD1i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
68 fFMD1i->SetXTitle("Hits (MC)");
69 fFMD2i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
70 fFMD2i->SetXTitle("Hits (MC)");
71 fFMD2o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
72 fFMD2o->SetXTitle("Hits (MC)");
73 fFMD3i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
74 fFMD3i->SetXTitle("Hits (MC)");
75 fFMD3o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
76 fFMD3o->SetXTitle("Hits (MC)");
77 fFMD1i->SetDirectory(0);
78 fFMD2i->SetDirectory(0);
79 fFMD2o->SetDirectory(0);
80 fFMD3i->SetDirectory(0);
81 fFMD3o->SetDirectory(0);
82 fSumEta = new TH1D("mcSumEta", "MC INEL Truth", 200, -4, 6);
83 fSumEta->SetXTitle("#eta");
84 fSumEta->SetYTitle("dN_{ch}/d#eta");
85 fSumEta->SetDirectory(0);
86 fSumEta->Sumw2();
87 fSumEta->SetMarkerColor(kOrange+2);
88 fSumEta->SetMarkerStyle(22);
89 fSumEta->SetFillColor(0);
90 fSumEta->SetFillStyle(0);
5bb5d1f6 91
92 fOper = new AliFMDFloatMap(0,0,0,0);
93 fOperComp = new TH2I("operComp", "Operation vs # track refs",
94 kMergedInto, kNone-.5, kMergedInto+.5,
95 20, -.5, 19.5);
96 fOperComp->SetXTitle("Operation");
97 fOperComp->SetYTitle("# of track refs in sector");
98 fOperComp->SetZTitle("Observations");
99 fOperComp->GetXaxis()->SetBinLabel(kNone, "None");
100 fOperComp->GetXaxis()->SetBinLabel(kCandidate, "Candidate");
101 fOperComp->GetXaxis()->SetBinLabel(kMergedWithOther, "Merged w/other");
102 fOperComp->GetXaxis()->SetBinLabel(kMergedInto, "Merged into");
103 fOperComp->SetDirectory(0);
dc0b1641 104
5bb5d1f6 105 fThetaVsNr = new TH2D("thetaVsNr", "#theta of track vs # track references",
106 360, 0, 360, 20, -.5, 19.5);
107 fThetaVsNr->SetXTitle("#theta [degrees]");
108 fThetaVsNr->SetYTitle("# of track references");
109 fThetaVsNr->SetDirectory(0);
dc0b1641 110}
111
112//____________________________________________________________________
113AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o)
114 : AliFMDSharingFilter(o),
115 fFMD1i(o.fFMD1i),
116 fFMD2i(o.fFMD2i),
117 fFMD2o(o.fFMD2o),
118 fFMD3i(o.fFMD3i),
119 fFMD3o(o.fFMD3o),
5bb5d1f6 120 fSumEta(o.fSumEta),
121 fOperComp(o.fOperComp),
122 fThetaVsNr(o.fThetaVsNr),
123 fOnlyPrimary(o.fOnlyPrimary)
dc0b1641 124{
7984e5f7 125 //
126 // Copy constructor
127 //
128 // Parameters:
129 // o Object to copy from
130 //
dc0b1641 131}
132
133//____________________________________________________________________
134AliFMDMCSharingFilter::~AliFMDMCSharingFilter()
135{
7984e5f7 136 //
137 // Destructor
138 //
dc0b1641 139 if (fFMD1i) delete fFMD1i;
140 if (fFMD2i) delete fFMD2i;
141 if (fFMD2o) delete fFMD2o;
142 if (fFMD3i) delete fFMD3i;
143 if (fFMD3o) delete fFMD3o;
144 if (fSumEta) delete fSumEta;
145}
146
147//____________________________________________________________________
148AliFMDMCSharingFilter&
149AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o)
150{
7984e5f7 151 //
152 // Assignment operator
153 //
154 // Parameters:
155 // o Object to assign from
156 //
157 // Return:
158 // Reference to this
159 //
dc0b1641 160 AliFMDSharingFilter::operator=(o);
5bb5d1f6 161 fOnlyPrimary = o.fOnlyPrimary;
dc0b1641 162 return *this;
163}
164
165//____________________________________________________________________
166void
167AliFMDMCSharingFilter::StoreParticle(UShort_t d,
168 Char_t r,
169 UShort_t s,
170 UShort_t t,
5bb5d1f6 171 UShort_t nR,
172 Double_t theta,
dc0b1641 173 AliESDFMD& output) const
174{
7984e5f7 175 //
176 // Store a particle hit in FMD<i>dr</i>[<i>s,t</i>] in @a output
177 //
178 // Parameters:
179 // d Detector
180 // r Ring
181 // s Sector
182 // t Strip
5bb5d1f6 183 // nR Number of references to this particle in this sector
7984e5f7 184 // output Output ESD object
185 //
dc0b1641 186 Double_t old = output.Multiplicity(d,r,s,t);
187 if (old == AliESDFMD::kInvalidMult) old = 0;
5bb5d1f6 188 if (fOper) fOperComp->Fill(fOper->operator()(d,r,s,t), nR);
189 if (theta < 0) theta += 2*TMath::Pi();
190 theta *= 180. / TMath::Pi();
191 fThetaVsNr->Fill(theta, nR);
dc0b1641 192 output.SetMultiplicity(d,r,s,t,old+1);
193}
194
195//____________________________________________________________________
196Bool_t
197AliFMDMCSharingFilter::FilterMC(const AliESDFMD& input,
198 const AliMCEvent& event,
199 Double_t vz,
4cbdf467 200 AliESDFMD& output,
201 TH2D* primary)
dc0b1641 202{
7984e5f7 203 //
4cbdf467 204 // Filter the input kinematics and track references, using
205 // some of the ESD information
7984e5f7 206 //
207 // Parameters:
4cbdf467 208 // input Input ESD event
209 // event Input MC event
210 // vz Vertex position
211 // output Output ESD-like object
212 // primary Per-event histogram of primaries
213 //
7984e5f7 214 // Return:
4cbdf467 215 // True on succes, false otherwise
7984e5f7 216 //
dc0b1641 217 output.Clear();
218
5bb5d1f6 219 // Increase event count - stored in
220 // underflow bin
221 fSumEta->AddBinContent(0);
222
dc0b1641 223 // Copy eta values to output
224 for (UShort_t ed = 1; ed <= 3; ed++) {
225 UShort_t nq = (ed == 1 ? 1 : 2);
226 for (UShort_t eq = 0; eq < nq; eq++) {
227 Char_t er = (eq == 0 ? 'I' : 'O');
228 UShort_t ns = (eq == 0 ? 20 : 40);
229 UShort_t nt = (eq == 0 ? 512 : 256);
230 for (UShort_t es = 0; es < ns; es++)
231 for (UShort_t et = 0; et < nt; et++)
232 output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
233 }
234 }
235 AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
e58000b7 236 Int_t nTracks = stack->GetNtrack();//event.GetNumberOfTracks();
237 Int_t nPrim = stack->GetNprimary();//event.GetNumberOfPrimary();
dc0b1641 238 for (Int_t iTr = 0; iTr < nTracks; iTr++) {
239 AliMCParticle* particle =
240 static_cast<AliMCParticle*>(event.GetTrack(iTr));
241
242 // Check the returned particle
243 if (!particle) continue;
244
245 // Check if this charged and a primary
246 Bool_t isCharged = particle->Charge() != 0;
247 if (!isCharged) continue;
248 Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
249
5bb5d1f6 250 // Pseudo rapidity and azimuthal angle
251 Double_t eta = particle->Eta();
252 Double_t phi = particle->Phi();
253
dc0b1641 254 // Fill 'dn/deta' histogram
e58000b7 255 if (isPrimary && iTr < nPrim) {
5bb5d1f6 256 // Avoid under count - used to store event count
257 if (eta >= fSumEta->GetXaxis()->GetXmin()) fSumEta->Fill(eta);
258 primary->Fill(eta, phi);
4cbdf467 259 }
dc0b1641 260
5bb5d1f6 261 // Bail out if we're only processing primaries - perhaps we should
262 // track back to the original primary?
263 if (fOnlyPrimary && !isPrimary) continue;
264
dc0b1641 265 Int_t nTrRef = particle->GetNumberOfTrackReferences();
266 Int_t longest = -1;
267 Double_t angle = 0;
268 UShort_t oD = 0, oS = 1024, oT = 1024;
269 Char_t oR = '\0';
5bb5d1f6 270 UShort_t nC = 0;
271 Double_t oTheta = 0;
dc0b1641 272 for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
273 AliTrackReference* ref = particle->GetTrackReference(iTrRef);
274
275 // Check existence
276 if (!ref) continue;
277
278 // Check that we hit an FMD element
279 if (ref->DetectorId() != AliTrackReference::kFMD)
280 continue;
281
5bb5d1f6 282 // Count number of track refs in this sector
283 nC++;
284
dc0b1641 285 // Get the detector coordinates
286 UShort_t d, s, t;
287 Char_t r;
288 AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
5bb5d1f6 289 // If this is a new detector/ring, then reset the other one
290 if (oD > 0 && oR != '\0' && oS != 1024 &&
291 (d != oD || r != oR || s != oS)) {
dc0b1641 292 longest = -1;
5bb5d1f6 293 angle = 0;
294 StoreParticle(oD, oR, oS, oT, nC, oTheta, output);
295 nC = 0;
296 oD = 0;
297 oR = '\0';
298 oS = 1024;
299 oT = 1024;
dc0b1641 300 }
301
302 oD = d;
303 oR = r;
304 oS = s;
305 oT = t;
306
307 // The longest passage is determined through the angle
308 Double_t x = ref->X();
309 Double_t y = ref->Y();
310 Double_t z = ref->Z()-vz;
311 Double_t rr = TMath::Sqrt(x*x+y*y);
312 Double_t theta= TMath::ATan2(rr,z);
313 Double_t ang = TMath::Abs(TMath::Pi()-theta);
314 if (ang > angle) {
315 longest = iTrRef;
316 angle = ang;
317 }
5bb5d1f6 318 oTheta = theta;
dc0b1641 319 } // Loop over track references
320 if (longest < 0) continue;
321
322 // Get the reference corresponding to the longest path through the detector
323 AliTrackReference* ref = particle->GetTrackReference(longest);
324
325 // Get the detector coordinates
326 UShort_t d, s, t;
327 Char_t r;
328 AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
329
5bb5d1f6 330 StoreParticle(d,r,s,t,nC,particle->Theta(),output);
dc0b1641 331 } // Loop over tracks
332 return kTRUE;
333}
334
335//____________________________________________________________________
336void
337AliFMDMCSharingFilter::CompareResults(const AliESDFMD& esd,
338 const AliESDFMD& mc)
339{
7984e5f7 340 //
341 // Compare the result of merging to the monte-carlo truth. This
342 // fills the correlation histograms
343 //
344 // Parameters:
345 // esd ESD after sharing correction
346 // mc MC ESD
347 //
348
dc0b1641 349 // Copy eta values to output
350 for (UShort_t d = 1; d <= 3; d++) {
351 UShort_t nq = (d == 1 ? 1 : 2);
352 for (UShort_t q = 0; q < nq; q++) {
353 Char_t r = (q == 0 ? 'I' : 'O');
354 UShort_t ns = (q == 0 ? 20 : 40);
355 UShort_t nt = (q == 0 ? 512 : 256);
356 TH2* co = 0;
357 switch (d) {
358 case 1: co = fFMD1i; break;
359 case 2: co = (q == 0 ? fFMD2i : fFMD2o); break;
360 case 3: co = (q == 0 ? fFMD3i : fFMD3o); break;
361 }
362
363 for (UShort_t s = 0; s < ns; s++) {
364 for (UShort_t t = 0; t < nt; t++) {
365 Float_t mEsd = esd.Multiplicity(d, r, s, t);
366 Float_t mMc = mc.Multiplicity(d, r, s, t);
367
368 co->Fill(mMc, mEsd);
369 }
370 }
371 }
372 }
373}
374
375//____________________________________________________________________
376void
377AliFMDMCSharingFilter::DefineOutput(TList* dir)
378{
7984e5f7 379 //
380 // Define the output histograms. These are put in a sub list of the
381 // passed list. The histograms are merged before the parent task calls
382 // AliAnalysisTaskSE::Terminate
383 //
384 // Parameters:
385 // dir Directory to add to
386 //
dc0b1641 387 AliFMDSharingFilter::DefineOutput(dir);
388 TList* d = static_cast<TList*>(dir->FindObject(GetName()));
389 TList* cd = new TList;
390 cd->SetName("esd_mc_comparion");
391 d->Add(cd);
392 cd->Add(fFMD1i);
393 cd->Add(fFMD2i);
394 cd->Add(fFMD2o);
395 cd->Add(fFMD3i);
396 cd->Add(fFMD3o);
397 dir->Add(fSumEta);
5bb5d1f6 398 cd->Add(fOperComp);
399 cd->Add(fThetaVsNr);
dc0b1641 400}
401
402//____________________________________________________________________
403void
fb3430ac 404AliFMDMCSharingFilter::ScaleHistograms(const TList* dir, Int_t nEvents)
dc0b1641 405{
7984e5f7 406 //
407 // Scale the histograms to the total number of events
408 //
409 // Parameters:
410 // dir Where the output is
411 // nEvents Number of events
412 //
dc0b1641 413 AliFMDSharingFilter::ScaleHistograms(dir, nEvents);
414 TH1D* sumEta = static_cast<TH1D*>(dir->FindObject("mcSumEta"));
415 if (!sumEta) {
416 AliWarning(Form("No mcSumEta histogram found in output list"));
417 return;
418 }
5bb5d1f6 419 Double_t n = nEvents; // sumEta->GetBinContent(0);
420 sumEta->Scale(1. / n, "width");
dc0b1641 421}
422
423//____________________________________________________________________
424void
425AliFMDMCSharingFilter::Print(Option_t* option) const
426{
7984e5f7 427 //
428 // Print information
429 //
430 // Parameters:
431 // option Not used
432 //
dc0b1641 433 char ind[gROOT->GetDirLevel()+1];
434 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
435 ind[gROOT->GetDirLevel()] = '\0';
436 AliFMDSharingFilter::Print(option);
5bb5d1f6 437 std::cout << std::boolalpha
438 << ind << " Only primary tracks: " << fOnlyPrimary
439 << std::noboolalpha << std::endl;
dc0b1641 440}
441
442//____________________________________________________________________
443//
444// EOF
445//