Refactored Ali{SPD,FMD}MCTrackDensity into a common base class.
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDMCTrackDensity.cxx
1 #include "AliFMDMCTrackDensity.h"
2 #include <AliESDFMD.h>
3 #include <AliTrackReference.h>
4 #include <TMath.h>
5 #include "AliFMDStripIndex.h"
6 #include <AliLog.h>
7 #include <TH2D.h>
8 #include <TH1D.h>
9 #include <TList.h>
10 #include <TROOT.h>
11 #include <iostream>
12
13 //____________________________________________________________________
14 void
15 AliFMDMCTrackDensity::State::Clear(Bool_t alsoCount)
16 {
17   angle       = 0;
18   oldDetector = 0;
19   oldRing     = '\0';
20   oldSector   = 1024;
21   oldStrip    = 1024;
22   startStrip  = 1024;
23   nRefs       = 0;
24   nStrips     = 0;
25   longest     = 0x0;
26   if (alsoCount) count = 0;
27 }
28
29 //____________________________________________________________________
30 AliFMDMCTrackDensity::State&
31 AliFMDMCTrackDensity::State::operator=(const State& o)
32 {
33   angle          = o.angle;
34   oldDetector    = o.oldDetector;
35   oldRing        = o.oldRing;
36   oldSector      = o.oldSector;
37   oldStrip       = o.oldStrip;
38   startStrip     = o.startStrip;
39   nRefs          = o.nRefs;
40   nStrips        = o.nStrips;
41   count          = o.count;
42   longest        = o.longest;
43   return *this;
44 }
45
46 //____________________________________________________________________
47 AliFMDMCTrackDensity::AliFMDMCTrackDensity()
48   : AliBaseMCTrackDensity(), 
49     fState(),
50     fMaxConsequtiveStrips(3), 
51     fNr(0), 
52     fNt(0), 
53     fNc(0),
54     fNcr(0),
55     fOutput(0)
56 {
57   // Default constructor 
58 }
59
60 //____________________________________________________________________
61 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
62   : AliBaseMCTrackDensity("fmdMCTrackDensity"), 
63     fState(),
64     fMaxConsequtiveStrips(3), 
65     fNr(0), 
66     fNt(0), 
67     fNc(0),
68     fNcr(0),
69     fOutput(0)
70 {
71   // Normal constructor constructor 
72 }
73
74 //____________________________________________________________________
75 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
76   : AliBaseMCTrackDensity(o),
77     fState(o.fState),
78     fMaxConsequtiveStrips(o.fMaxConsequtiveStrips), 
79     fNr(o.fNr), 
80     fNt(o.fNt), 
81     fNc(o.fNc),
82     fNcr(o.fNcr),
83     fOutput(o.fOutput)
84 {
85   // Normal constructor constructor 
86 }
87
88 //____________________________________________________________________
89 AliFMDMCTrackDensity&
90 AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
91 {
92   // Assignment operator 
93   if (&o == this) return *this; 
94   AliBaseMCTrackDensity::operator=(o);
95   fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
96   fNr                   = o.fNr;
97   fNt                   = o.fNt;
98   fNc                   = o.fNc;
99   fNcr                  = o.fNcr;
100   fState                = o.fState;
101   fOutput               = o.fOutput;
102
103   return *this;
104 }
105
106 //____________________________________________________________________
107 void
108 AliFMDMCTrackDensity::DefineOutput(TList* l)
109 {
110   AliBaseMCTrackDensity::DefineOutput(l);
111   TList* ll = static_cast<TList*>(l->FindObject(GetTitle()));
112   if (!ll) ll = l;
113
114   fNr = new TH1D("clusterRefs", "# track references per cluster",
115                  21, -.5, 20.5);
116   fNr->SetXTitle("# of track references/cluster");
117   fNr->SetFillColor(kRed+1);
118   fNr->SetFillStyle(3001);
119   fNr->SetDirectory(0);
120   ll->Add(fNr);
121
122   fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
123   fNt->SetXTitle("Cluster size [strips]");
124   fNt->SetFillColor(kBlue+1);
125   fNt->SetFillStyle(3001);
126   fNt->SetDirectory(0);
127   ll->Add(fNt);
128
129   fNc = new TH1D("nClusters", "# clusters per track", 21, -.5, 20.5);
130   fNc->SetXTitle("# clusters");
131   fNc->SetFillColor(kGreen+1);
132   fNc->SetFillStyle(3001);
133   fNc->SetDirectory(0);
134   ll->Add(fNc);
135
136   fNcr = new TH2D("clusterVsRefs", "# clusters vs. # refs", 
137                   21, -.5, 20.5, 21, -.5, 20.5);
138   fNcr->SetXTitle("# References");
139   fNcr->SetYTitle("# Clusters");
140   fNcr->SetOption("COLZ");
141   fNcr->SetDirectory(0);
142   ll->Add(fNcr);
143   
144                   
145 }
146 //____________________________________________________________________
147 Int_t
148 AliFMDMCTrackDensity::GetDetectorId() const
149 {
150   return AliTrackReference::kFMD;
151 }
152
153 //____________________________________________________________________
154 void
155 AliFMDMCTrackDensity::BeginTrackRefs()
156 {
157   fState.Clear(true);
158 }
159
160 //____________________________________________________________________
161 void
162 AliFMDMCTrackDensity::EndTrackRefs(Int_t nRefs)
163 {
164   fNc->Fill(fState.count);
165   fNcr->Fill(nRefs, fState.count);
166   fState.Clear(true);
167 }
168   
169 //____________________________________________________________________
170 AliTrackReference*
171 AliFMDMCTrackDensity::ProcessRef(AliMCParticle*       particle,
172                                  const AliMCParticle* mother,
173                                  AliTrackReference*   ref)
174 {
175   // Get the detector coordinates 
176   UShort_t d, s, t;
177   Char_t r;
178   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
179     
180   // Calculate distance of previous reference to base of cluster 
181   UShort_t nT = TMath::Abs(t - fState.startStrip) + 1;
182
183   // Now check if we should flush to output 
184   Bool_t used = false;
185
186   // If this is a new detector/ring, then reset the other one 
187   // Check if we have a valid old detectorm ring, and sector 
188   if (fState.oldDetector >  0 && 
189       fState.oldRing     != '\0' && 
190       fState.oldSector   != 1024) {
191     // New detector, new ring, or new sector 
192     if (d != fState.oldDetector   || 
193         r != fState.oldRing       || 
194         s != fState.oldSector) {
195       if (fDebug) Info("Process", "New because new sector");
196       used = true;
197     }
198     else if (nT > fMaxConsequtiveStrips) {
199       if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)", 
200                        fState.nStrips, t, fState.oldStrip, fState.startStrip);
201       used = true;
202     }
203   }
204   if (used) {
205     if (fDebug) 
206       Info("Process", "I=%p L=%p D=%d (was %d), R=%c (was %c), "
207            "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
208            ref, fState.longest, 
209            d, fState.oldDetector, 
210            r, fState.oldRing, 
211            s, fState.oldSector, 
212            t, fState.oldStrip, 
213            fState.nStrips, fMaxConsequtiveStrips);
214     // Int_t nnT   = TMath::Abs(fState.oldStrip - fState.startStrip) + 1;
215     StoreParticle(particle, mother, fState.longest);
216     fState.Clear(false);
217   }
218
219   // If base of cluster not set, set it here. 
220   if (fState.startStrip == 1024) fState.startStrip = t;
221   
222   // Calculate distance of previous reference to base of cluster 
223   fState.nStrips = TMath::Abs(t - fState.startStrip) + 1;
224
225   // Count number of track refs in this sector 
226   fState.nRefs++;
227
228   fState.oldDetector = d;
229   fState.oldRing     = r;
230   fState.oldSector   = s;
231   fState.oldStrip    = t;
232
233   // Debug output 
234   if (fDebug) {
235     if (t == fState.startStrip) 
236       Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]", 
237            d, r, s, t);
238     else 
239       Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
240            "length=%3d (now in %3d, previous %3d)", 
241            d, r, s, fState.startStrip, fState.nStrips, t, fState.oldStrip);
242   }
243     
244   // The longest passage is determined through the angle 
245   Double_t ang  = GetTrackRefTheta(ref);
246   if (ang > fState.angle) {
247     fState.longest = ref;
248     fState.angle   = ang;
249   }
250
251   return fState.longest;
252 }
253
254 //____________________________________________________________________
255 Double_t
256 AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle, 
257                                     const AliMCParticle* mother, 
258                                     AliTrackReference*   ref) const
259 {
260   Double_t w = 
261     AliBaseMCTrackDensity::StoreParticle(particle, mother, ref);
262   if (w <= 0) return w;
263
264   // Get the detector coordinates 
265   UShort_t d, s, t;
266   Char_t r;
267   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
268
269   // Check if we have value already 
270   Double_t old = fOutput->Multiplicity(d,r,s,t);
271
272   // If invalid, force it valid 
273   if (old == AliESDFMD::kInvalidMult) old = 0;
274
275   // Increment count 
276   fOutput->SetMultiplicity(d,r,s,t,old+w);
277
278   // Fill histograms 
279   fNr->Fill(fState.nRefs);
280   fNt->Fill(fState.nStrips);
281
282   fState.count++;
283
284   return w;
285 }  
286
287 //____________________________________________________________________
288 Bool_t
289 AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input, 
290                                 const AliMCEvent& event,
291                                 Double_t          vz,
292                                 AliESDFMD&        output, 
293                                 TH2D*             primary)
294 {
295   // 
296   // Filter the input kinematics and track references, using 
297   // some of the ESD information
298   // 
299   // Parameters:
300   //    input   Input ESD event
301   //    event   Input MC event
302   //    vz      Vertex position 
303   //    output  Output ESD-like object
304   //    primary Per-event histogram of primaries 
305   //
306   // Return:
307   //    True on succes, false otherwise 
308   //
309   output.Clear();
310   fOutput = &output;
311
312   // Copy eta values to output 
313   for (UShort_t ed = 1; ed <= 3; ed++) { 
314     UShort_t nq = (ed == 1 ? 1 : 2);
315     for (UShort_t eq = 0; eq < nq; eq++) {
316       Char_t   er = (eq == 0 ? 'I' : 'O');
317       UShort_t ns = (eq == 0 ?  20 :  40);
318       UShort_t nt = (eq == 0 ? 512 : 256);
319       for (UShort_t es = 0; es < ns; es++) 
320         for (UShort_t et = 0; et < nt; et++) 
321           output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
322     }
323   }
324
325   return ProcessTracks(event, vz, primary);
326 }
327 //____________________________________________________________________
328 void
329 AliFMDMCTrackDensity::Print(Option_t* option) const 
330 {
331   AliBaseMCTrackDensity::Print(option);
332   char ind[gROOT->GetDirLevel()+1];
333   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
334   ind[gROOT->GetDirLevel()] = '\0';
335   std::cout << ind << " Max cluster size:       " << fMaxConsequtiveStrips
336             << std::endl;
337   
338 }
339
340 //____________________________________________________________________
341 //
342 // EOF
343 //