Added sub-tasks to do hit counting based on track-refs.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDMCTrackDensity.cxx
CommitLineData
f68d9069 1#include "AliFMDMCTrackDensity.h"
2#include <AliESDFMD.h>
3#include <AliMCEvent.h>
4#include <AliTrackReference.h>
5#include <AliStack.h>
6#include <TMath.h>
7#include "AliFMDStripIndex.h"
8#include "AliFMDFloatMap.h"
9#include <AliLog.h>
10#include <TH2D.h>
11#include <TH1D.h>
12#include <TList.h>
13#include <TROOT.h>
14#include <iostream>
15
16//____________________________________________________________________
17AliFMDMCTrackDensity::AliFMDMCTrackDensity()
18 : TNamed(),
19 fUseOnlyPrimary(false),
20 fMaxConsequtiveStrips(3),
21 fNr(0),
22 fNt(0),
23 fBinFlow(0),
24 fEtaBinFlow(0),
25 fPhiBinFlow(0),
26 fDebug(false)
27{
28 // Default constructor
29}
30
31//____________________________________________________________________
32AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
33 : TNamed("fmdMCTrackDensity","fmdMCTrackDensity"),
34 fUseOnlyPrimary(false),
35 fMaxConsequtiveStrips(3),
36 fNr(0),
37 fNt(0),
38 fBinFlow(0),
39 fEtaBinFlow(0),
40 fPhiBinFlow(0),
41 fDebug(false)
42{
43 // Normal constructor constructor
44}
45
46//____________________________________________________________________
47AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
48 : TNamed(o),
49 fUseOnlyPrimary(o.fUseOnlyPrimary),
50 fMaxConsequtiveStrips(o.fMaxConsequtiveStrips),
51 fNr(o.fNr),
52 fNt(o.fNt),
53 fBinFlow(o.fBinFlow),
54 fEtaBinFlow(o.fEtaBinFlow),
55 fPhiBinFlow(o.fPhiBinFlow),
56 fDebug(o.fDebug)
57{
58 // Normal constructor constructor
59}
60
61//____________________________________________________________________
62AliFMDMCTrackDensity&
63AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
64{
65 // Assignment operator
66 TNamed::operator=(o);
67 fUseOnlyPrimary = o.fUseOnlyPrimary;
68 fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
69 fNr = o.fNr;
70 fNt = o.fNt;
71 fBinFlow = o.fBinFlow;
72 fEtaBinFlow = o.fEtaBinFlow;
73 fPhiBinFlow = o.fPhiBinFlow;
74 fDebug = o.fDebug;
75 return *this;
76}
77
78//____________________________________________________________________
79void
80AliFMDMCTrackDensity::DefineOutput(TList* l)
81{
82 fNr = new TH1D("clusterRefs", "# track references per cluster",
83 21, -.5, 20.5);
84 fNr->SetXTitle("# of track references/cluster");
85 fNr->SetDirectory(0);
86 l->Add(fNr);
87
88 fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
89 fNt->SetXTitle("Cluster size [strips]");
90 fNt->SetDirectory(0);
91 l->Add(fNt);
92
93 fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
94 200, -5, 5, 40, -180, 180);
95 fBinFlow->SetXTitle("#Delta#eta");
96 fBinFlow->SetYTitle("#Delta#varphi");
97 fBinFlow->SetDirectory(0);
98 l->Add(fBinFlow);
99
100 fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta",
101 200, -4, 6, 200, -5, 5);
102 fEtaBinFlow->SetXTitle("#eta of strip");
103 fEtaBinFlow->SetYTitle("#Delta#eta");
104 fEtaBinFlow->SetDirectory(0);
105 l->Add(fEtaBinFlow);
106
107 fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi",
108 40, 0, 360, 40, -180, 180);
109 fPhiBinFlow->SetXTitle("#varphi of strip");
110 fPhiBinFlow->SetYTitle("#Delta#varphi");
111 fPhiBinFlow->SetDirectory(0);
112 l->Add(fPhiBinFlow);
113}
114
115//____________________________________________________________________
116void
117AliFMDMCTrackDensity::StoreParticle(AliMCParticle* particle,
118 const AliMCParticle* mother,
119 Int_t longest,
120 Double_t vz,
121 UShort_t nC,
122 UShort_t nT,
123 AliESDFMD& output) const
124{
125 // Store a particle.
126 if (longest < 0) return;
127
128 AliTrackReference* ref = particle->GetTrackReference(longest);
129 if (!ref) return;
130
131 // Get the detector coordinates
132 UShort_t d, s, t;
133 Char_t r;
134 AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
135
136 // Check if we have value already
137 Double_t old = output.Multiplicity(d,r,s,t);
138
139 // If invalid, force it valid
140 if (old == AliESDFMD::kInvalidMult) old = 0;
141
142 // Increment count
143 output.SetMultiplicity(d,r,s,t,old+1);
144
145 // Get track-reference stuff
146 Double_t x = ref->X();
147 Double_t y = ref->Y();
148 Double_t z = ref->Z()-vz;
149 Double_t rr = TMath::Sqrt(x*x+y*y);
150 Double_t thetaR = TMath::ATan2(rr,z);
151 Double_t phiR = TMath::ATan2(y,x);
152 Double_t etaR = -TMath::Log(TMath::Tan(thetaR/2));
153
154 // Correct angle and convert to degrees
155 if (thetaR < 0) thetaR += 2*TMath::Pi();
156 thetaR *= 180. / TMath::Pi();
157 if (phiR < 0) phiR += 2*TMath::Pi();
158 phiR *= 180. / TMath::Pi();
159
160 // Fill histograms
161 fNr->Fill(nC);
162 fNt->Fill(nT);
163
164 const AliMCParticle* mp = (mother ? mother : particle);
165 Double_t dEta = mp->Eta() - etaR;
166 Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
167 if (dPhi > 180) dPhi -= 360;
168 if (dPhi < -180) dPhi += 360;
169 fBinFlow->Fill(dEta, dPhi);
170 fEtaBinFlow->Fill(etaR, dEta);
171 fPhiBinFlow->Fill(phiR, dPhi);
172
173 if (fDebug)
174 Info("StoreParticle", "Store @ FMD%d%c[%2d,%3d] nStrips=%3d, "
175 "dEta=%7.4f, dPhi=%d",
176 d, r, s, t, nT, dEta, int(dPhi+.5));
177}
178
179//____________________________________________________________________
180Double_t
181AliFMDMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref,
182 Double_t vz) const
183{
184 // Get the incidient angle of the track reference.
185 Double_t x = ref->X();
186 Double_t y = ref->Y();
187 Double_t z = ref->Z()-vz;
188 Double_t rr = TMath::Sqrt(x*x+y*y);
189 Double_t theta= TMath::ATan2(rr,z);
190 Double_t ang = TMath::Abs(TMath::Pi()-theta);
191 return ang;
192}
193
194//____________________________________________________________________
195const AliMCParticle*
196AliFMDMCTrackDensity::GetMother(Int_t iTr,
197 const AliMCEvent& event) const
198{
199 //
200 // Track down primary mother
201 //
202 Int_t i = iTr;
203 do {
204 const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
205 if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
206
207 i = p->GetMother();
208 } while (i > 0);
209
210 return 0;
211}
212
213//____________________________________________________________________
214Bool_t
215AliFMDMCTrackDensity::Calculate(const AliESDFMD& input,
216 const AliMCEvent& event,
217 Double_t vz,
218 AliESDFMD& output,
219 TH2D* primary)
220{
221 //
222 // Filter the input kinematics and track references, using
223 // some of the ESD information
224 //
225 // Parameters:
226 // input Input ESD event
227 // event Input MC event
228 // vz Vertex position
229 // output Output ESD-like object
230 // primary Per-event histogram of primaries
231 //
232 // Return:
233 // True on succes, false otherwise
234 //
235 output.Clear();
236
237 // Copy eta values to output
238 for (UShort_t ed = 1; ed <= 3; ed++) {
239 UShort_t nq = (ed == 1 ? 1 : 2);
240 for (UShort_t eq = 0; eq < nq; eq++) {
241 Char_t er = (eq == 0 ? 'I' : 'O');
242 UShort_t ns = (eq == 0 ? 20 : 40);
243 UShort_t nt = (eq == 0 ? 512 : 256);
244 for (UShort_t es = 0; es < ns; es++)
245 for (UShort_t et = 0; et < nt; et++)
246 output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
247 }
248 }
249 AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
250 Int_t nTracks = stack->GetNtrack();//event.GetNumberOfTracks();
251 Int_t nPrim = stack->GetNprimary();//event.GetNumberOfPrimary();
252 for (Int_t iTr = 0; iTr < nTracks; iTr++) {
253 AliMCParticle* particle =
254 static_cast<AliMCParticle*>(event.GetTrack(iTr));
255
256 // Check the returned particle
257 if (!particle) continue;
258
259 // Check if this charged and a primary
260 Bool_t isCharged = particle->Charge() != 0;
261 if (!isCharged) continue;
262 Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
263
264 // Pseudo rapidity and azimuthal angle
265 Double_t eta = particle->Eta();
266 Double_t phi = particle->Phi();
267
268 // Fill 'dn/deta' histogram
269 if (isPrimary && iTr < nPrim) {
270 if (primary) primary->Fill(eta, phi);
271 }
272
273 // Bail out if we're only processing primaries - perhaps we should
274 // track back to the original primary?
275 if (fUseOnlyPrimary && !isPrimary) continue;
276
277 Int_t nTrRef = particle->GetNumberOfTrackReferences();
278 Int_t longest = -1;
279 Double_t angle = 0;
280 UShort_t oD = 0, oS = 1024, oT = 1024, ooT=1024;
281 Char_t oR = '\0';
282 UShort_t nC = 0;
283 UShort_t nT = 0;
284 Double_t oTheta = 0;
285 for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
286 AliTrackReference* ref = particle->GetTrackReference(iTrRef);
287
288 // Check existence
289 if (!ref) continue;
290
291 // Check that we hit an FMD element
292 if (ref->DetectorId() != AliTrackReference::kFMD)
293 continue;
294
295 // Get the detector coordinates
296 UShort_t d, s, t;
297 Char_t r;
298 AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
299
300 // Calculate length of last and second to last strips.
301
302 // IF base of cluster not set, set it here.
303 if (ooT == 1024) ooT = t;
304
305 // Calculate distance of previous reference to base of cluster
306 nT = TMath::Abs(t - ooT) + 1;
307
308 // Count number of track refs in this sector
309 nC++;
310
311 Bool_t used = false;
312 // If this is a new detector/ring, then reset the other one
313 // Check if we have a valid old detectorm ring, and sector
314 if (oD > 0 && oR != '\0' && oS != 1024) {
315 // New detector, new ring, or new sector
316 if (d != oD || r != oR || s != oS) {
317 if (fDebug) Info("Process", "New because new sector");
318 used = true;
319 }
320 }
321 if (nT > fMaxConsequtiveStrips) {
322 if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
323 nT, t, oT, ooT);
324 used = true;
325 }
326 if (used) {
327 if (fDebug)
328 Info("Process", "I=%3d L=%3d D=%d (was %d), R=%c (was %c), "
329 "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
330 iTr, longest, d, oD, r, oR, s, oS, t, oT, nT,
331 fMaxConsequtiveStrips);
332 Int_t nnT = TMath::Abs(oT - ooT) + 1;
333 const AliMCParticle* mother = GetMother(iTr, event);
334 StoreParticle(particle, mother, longest, vz, nC, nnT, output);
335 longest = -1;
336 angle = 0;
337 nC = 1; // Reset track-ref counter - we have this->1
338 nT = 0; // Reset cluster size
339 oD = 0; // Reset detector
340 oR = '\0'; // Reset ring
341 oS = 1024; // Reset sector
342 oT = 1024; // Reset old strip
343 ooT = t; // Reset base
344 }
345 else if (fDebug) {
346 if (ooT == t)
347 Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
348 d, r, s, t);
349 else
350 Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
351 "length=%3d (now in %3d, previous %3d)",
352 d, r, s, ooT, nT, t, oT);
353 }
354
355 oD = d;
356 oR = r;
357 oS = s;
358 oT = t;
359 nT = TMath::Abs(t-ooT)+1;
360
361 // The longest passage is determined through the angle
362 Double_t ang = GetTrackRefTheta(ref, vz);
363 if (ang > angle) {
364 longest = iTrRef;
365 angle = ang;
366 }
367 oTheta = ang;
368 } // Loop over track references
369 if (longest < 0) continue; // Nothing found
370
371 // Get the reference corresponding to the longest path through the detector
372 if (fDebug)
373 Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)",
374 iTr, longest, nT, fMaxConsequtiveStrips);
375 const AliMCParticle* mother = GetMother(iTr, event);
376 StoreParticle(particle, mother, longest, vz, nC, nT, output);
377 } // Loop over tracks
378 return kTRUE;
379}
380//____________________________________________________________________
381void
382AliFMDMCTrackDensity::Print(Option_t* /*option*/) const
383{
384 char ind[gROOT->GetDirLevel()+1];
385 for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
386 ind[gROOT->GetDirLevel()] = '\0';
387 std::cout << ind << ClassName() << ": " << GetName() << '\n'
388 << std::boolalpha
389 << ind << " Only primary tracks: " << fUseOnlyPrimary << '\n'
390 << ind << " Max cluster size: " << fMaxConsequtiveStrips
391 << std::noboolalpha << std::endl;
392
393}
394
395//____________________________________________________________________
396//
397// EOF
398//