Fixed warnings [-Wunused-but-set-variable] from GCC 4.6 -
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDMCTrackDensity.cxx
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 //____________________________________________________________________
17 AliFMDMCTrackDensity::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 //____________________________________________________________________
32 AliFMDMCTrackDensity::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 //____________________________________________________________________
47 AliFMDMCTrackDensity::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 //____________________________________________________________________
62 AliFMDMCTrackDensity&
63 AliFMDMCTrackDensity::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 //____________________________________________________________________
79 void
80 AliFMDMCTrackDensity::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 //____________________________________________________________________
116 void
117 AliFMDMCTrackDensity::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 //____________________________________________________________________
180 Double_t
181 AliFMDMCTrackDensity::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 //____________________________________________________________________
195 const AliMCParticle*
196 AliFMDMCTrackDensity::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 //____________________________________________________________________
214 Bool_t
215 AliFMDMCTrackDensity::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 //____________________________________________________________________
381 void
382 AliFMDMCTrackDensity::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 //