8972f271277e5ac50f0ca3bdc50b0b72cd13a049
[u/mrichter/AliRoot.git] / PWGLF / 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 #include "AliGenHijingEventHeader.h"
16 #include <TF1.h>
17 #include <TGraph.h>
18
19 //____________________________________________________________________
20 AliFMDMCTrackDensity::AliFMDMCTrackDensity()
21   : TNamed(), 
22     fUseOnlyPrimary(false), 
23     fMaxConsequtiveStrips(3), 
24     fNr(0), 
25     fNt(0),
26     fBinFlow(0), 
27     fEtaBinFlow(0),
28     fPhiBinFlow(0),
29     fDebug(false)
30 {
31   // Default constructor 
32 }
33
34 //____________________________________________________________________
35 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
36   : TNamed("fmdMCTrackDensity","fmdMCTrackDensity"), 
37     fUseOnlyPrimary(false), 
38     fMaxConsequtiveStrips(3), 
39     fNr(0), 
40     fNt(0),
41     fBinFlow(0), 
42     fEtaBinFlow(0),
43     fPhiBinFlow(0),
44     fDebug(false)
45 {
46   // Normal constructor constructor 
47 }
48
49 //____________________________________________________________________
50 AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
51   : TNamed(o),
52     fUseOnlyPrimary(o.fUseOnlyPrimary), 
53     fMaxConsequtiveStrips(o.fMaxConsequtiveStrips), 
54     fNr(o.fNr), 
55     fNt(o.fNt),
56     fBinFlow(o.fBinFlow), 
57     fEtaBinFlow(o.fEtaBinFlow),
58     fPhiBinFlow(o.fPhiBinFlow),
59     fDebug(o.fDebug)
60 {
61   // Normal constructor constructor 
62 }
63
64 //____________________________________________________________________
65 AliFMDMCTrackDensity&
66 AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
67 {
68   // Assignment operator 
69   if (&o == this) return *this; 
70   TNamed::operator=(o);
71   fUseOnlyPrimary       = o.fUseOnlyPrimary;
72   fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
73   fNr                   = o.fNr;
74   fNt                   = o.fNt;
75   fBinFlow              = o.fBinFlow;
76   fEtaBinFlow           = o.fEtaBinFlow;
77   fPhiBinFlow           = o.fPhiBinFlow;
78   fDebug                = o.fDebug;
79   return *this;
80 }
81
82 //____________________________________________________________________
83 void
84 AliFMDMCTrackDensity::DefineOutput(TList* l)
85 {
86   fNr = new TH1D("clusterRefs", "# track references per cluster",
87                  21, -.5, 20.5);
88   fNr->SetXTitle("# of track references/cluster");
89   fNr->SetDirectory(0);
90   l->Add(fNr);
91
92   fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
93   fNt->SetXTitle("Cluster size [strips]");
94   fNt->SetDirectory(0);
95   l->Add(fNt);
96
97   fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow", 
98                       200, -5, 5, 40, -180, 180);
99   fBinFlow->SetXTitle("#Delta#eta");
100   fBinFlow->SetYTitle("#Delta#varphi");
101   fBinFlow->SetDirectory(0);
102   l->Add(fBinFlow);
103
104   fEtaBinFlow = new TH2D("binFlowEta", "#eta bin flow vs #eta", 
105                          200, -4, 6, 200, -5, 5);
106   fEtaBinFlow->SetXTitle("#eta of strip");
107   fEtaBinFlow->SetYTitle("#Delta#eta");
108   fEtaBinFlow->SetDirectory(0);
109   l->Add(fEtaBinFlow);
110
111   fPhiBinFlow = new TH2D("binFlowPhi", "#phi bin flow vs #phi", 
112                          40, 0, 360, 40, -180, 180);
113   fPhiBinFlow->SetXTitle("#varphi of strip");
114   fPhiBinFlow->SetYTitle("#Delta#varphi");
115   fPhiBinFlow->SetDirectory(0);
116   l->Add(fPhiBinFlow);
117 }
118
119 //____________________________________________________________________
120 void
121 AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle, 
122                                     const AliMCParticle* mother, 
123                                     Int_t                longest,
124                                     Double_t             vz,
125                                     UShort_t             nC, 
126                                     UShort_t             nT,
127                                     AliESDFMD&           output,
128                                     Double_t             w) const
129 {
130   // Store a particle. 
131   if (longest < 0) return;
132
133   AliTrackReference* ref = particle->GetTrackReference(longest);
134   if (!ref) return;
135     
136   // Get the detector coordinates 
137   UShort_t d, s, t;
138   Char_t r;
139   AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
140
141   // Check if we have value already 
142   Double_t old = output.Multiplicity(d,r,s,t);
143
144   // If invalid, force it valid 
145   if (old == AliESDFMD::kInvalidMult) old = 0;
146
147   // Increment count 
148   output.SetMultiplicity(d,r,s,t,old+w);
149
150   // Get track-reference stuff 
151   Double_t x      = ref->X();
152   Double_t y      = ref->Y();
153   Double_t z      = ref->Z()-vz;
154   Double_t rr     = TMath::Sqrt(x*x+y*y);
155   Double_t thetaR = TMath::ATan2(rr,z);
156   Double_t phiR   = TMath::ATan2(y,x);
157   Double_t etaR   = -TMath::Log(TMath::Tan(thetaR/2));
158   
159   // Correct angle and convert to degrees 
160   if (thetaR < 0) thetaR += 2*TMath::Pi();
161   thetaR *= 180. / TMath::Pi();
162   if (phiR < 0) phiR += 2*TMath::Pi();
163   phiR *= 180. / TMath::Pi();
164
165   // Fill histograms 
166   fNr->Fill(nC);
167   fNt->Fill(nT);
168
169   const AliMCParticle* mp = (mother ? mother : particle);
170   Double_t dEta = mp->Eta() - etaR;
171   Double_t dPhi = mp->Phi() * 180 / TMath::Pi() - phiR;
172   if (dPhi >  180) dPhi -= 360;
173   if (dPhi < -180) dPhi += 360;
174   fBinFlow->Fill(dEta, dPhi);
175   fEtaBinFlow->Fill(etaR, dEta);
176   fPhiBinFlow->Fill(phiR, dPhi);
177
178   if (fDebug) 
179     Info("StoreParticle", "Store @ FMD%d%c[%2d,%3d] nStrips=%3d, "
180          "dEta=%7.4f, dPhi=%d",
181          d, r, s, t, nT, dEta, int(dPhi+.5));
182 }
183
184 //____________________________________________________________________
185 Double_t
186 AliFMDMCTrackDensity::GetTrackRefTheta(const AliTrackReference* ref,
187                                        Double_t vz) const
188 {
189   // Get the incidient angle of the track reference. 
190   Double_t x    = ref->X();
191   Double_t y    = ref->Y();
192   Double_t z    = ref->Z()-vz;
193   Double_t rr   = TMath::Sqrt(x*x+y*y);
194   Double_t theta= TMath::ATan2(rr,z);
195   Double_t ang  = TMath::Abs(TMath::Pi()-theta);
196   return ang;
197 }
198                                     
199 //____________________________________________________________________
200 const AliMCParticle*
201 AliFMDMCTrackDensity::GetMother(Int_t     iTr,
202                                 const AliMCEvent& event) const
203 {
204   // 
205   // Track down primary mother 
206   // 
207   Int_t i  = iTr;
208   do { 
209     const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
210     if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
211     
212     i = p->GetMother();
213   } while (i > 0);
214
215   return 0;
216 }  
217
218 // #define USE_FLOW_WEIGHTS 1
219 //____________________________________________________________________
220 Bool_t
221 AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input, 
222                                 const AliMCEvent& event,
223                                 Double_t          vz,
224                                 AliESDFMD&        output, 
225                                 TH2D*             primary)
226 {
227   // 
228   // Filter the input kinematics and track references, using 
229   // some of the ESD information
230   // 
231   // Parameters:
232   //    input   Input ESD event
233   //    event   Input MC event
234   //    vz      Vertex position 
235   //    output  Output ESD-like object
236   //    primary Per-event histogram of primaries 
237   //
238   // Return:
239   //    True on succes, false otherwise 
240   //
241   output.Clear();
242
243   // Copy eta values to output 
244   for (UShort_t ed = 1; ed <= 3; ed++) { 
245     UShort_t nq = (ed == 1 ? 1 : 2);
246     for (UShort_t eq = 0; eq < nq; eq++) {
247       Char_t   er = (eq == 0 ? 'I' : 'O');
248       UShort_t ns = (eq == 0 ?  20 :  40);
249       UShort_t nt = (eq == 0 ? 512 : 256);
250       for (UShort_t es = 0; es < ns; es++) 
251         for (UShort_t et = 0; et < nt; et++) 
252           output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
253     }
254   }
255
256 #ifdef USE_FLOW_WEIGHTS
257   AliGenHijingEventHeader* hd = dynamic_cast<AliGenHijingEventHeader*>
258                 (event.GenEventHeader());
259   Double_t rp = (hd ? hd->ReactionPlaneAngle() : 0.);
260   Double_t b = (hd ? hd->ImpactParameter() : -1 );
261 #endif
262
263   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
264   Int_t nTracks   = stack->GetNtrack();//event.GetNumberOfTracks();
265   Int_t nPrim     = stack->GetNprimary();//event.GetNumberOfPrimary();
266   for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
267     AliMCParticle* particle = 
268       static_cast<AliMCParticle*>(event.GetTrack(iTr));
269     
270     // Check the returned particle 
271     if (!particle) continue;
272     
273     // Check if this charged and a primary 
274     Bool_t isCharged = particle->Charge() != 0;
275     if (!isCharged) continue;
276     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
277
278     // Pseudo rapidity and azimuthal angle 
279     Double_t eta = particle->Eta();
280     Double_t phi = particle->Phi();
281     
282     // Fill 'dn/deta' histogram 
283     if (isPrimary && iTr < nPrim) {
284       if (primary) primary->Fill(eta, phi);
285     }
286
287     // Bail out if we're only processing primaries - perhaps we should
288     // track back to the original primary?
289     if (fUseOnlyPrimary && !isPrimary) continue;
290
291     Int_t    nTrRef   = particle->GetNumberOfTrackReferences();
292     Int_t    longest  = -1;
293     Double_t angle    = 0;
294     UShort_t oD       = 0, oS = 1024, oT = 1024, ooT=1024;
295     Char_t   oR       = '\0';
296     UShort_t nC       = 0;
297     UShort_t nT       = 0;
298     // Double_t oTheta= 0;
299     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
300       AliTrackReference* ref = particle->GetTrackReference(iTrRef);
301       
302       // Check existence 
303       if (!ref) continue;
304
305       // Check that we hit an FMD element 
306       if (ref->DetectorId() != AliTrackReference::kFMD) 
307         continue;
308
309       // Get the detector coordinates 
310       UShort_t d, s, t;
311       Char_t r;
312       AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
313
314       // Calculate length of last and second to last strips. 
315
316       // IF base of cluster not set, set it here. 
317       if (ooT == 1024) ooT = t;
318
319       // Calculate distance of previous reference to base of cluster 
320       nT = TMath::Abs(t - ooT) + 1;
321
322       // Count number of track refs in this sector 
323       nC++;
324
325       Bool_t used = false;
326       // If this is a new detector/ring, then reset the other one 
327       // Check if we have a valid old detectorm ring, and sector 
328       if (oD >  0 && oR != '\0' && oS != 1024) {
329         // New detector, new ring, or new sector 
330         if (d != oD   || r != oR   || s != oS) {
331           if (fDebug) Info("Process", "New because new sector");
332           used = true;
333         }
334       }
335       if (nT > fMaxConsequtiveStrips) {
336         if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)", 
337                          nT, t, oT, ooT);
338         used = true;
339       }
340       if (used) {
341         if (fDebug) 
342           Info("Process", "I=%3d L=%3d D=%d (was %d), R=%c (was %c), "
343                "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
344                iTr, longest, d, oD, r, oR, s, oS, t, oT, nT, 
345                fMaxConsequtiveStrips);
346         Int_t nnT   = TMath::Abs(oT - ooT) + 1;
347         const AliMCParticle* mother = GetMother(iTr, event);
348         
349 #ifdef USE_FLOW_WEIGHTS
350         phi = (mother ? mother->Phi() : particle->Phi());
351         eta = (mother ? mother->Eta() : particle->Eta());
352         Double_t pt = (mother ? mother->Pt() : particle->Pt());
353         Int_t   id = (mother ? mother->PdgCode() : 2212);
354         Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
355 #else
356         Double_t weight = 1; // // 
357 #endif
358         StoreParticle(particle, mother, longest, vz, nC, nnT, output, weight);
359         longest = -1;
360         angle   = 0;
361         nC  = 1;    // Reset track-ref counter - we have this->1
362         nT  = 0;    // Reset cluster size 
363         oD  = 0;    // Reset detector
364         oR  = '\0'; // Reset ring 
365         oS  = 1024; // Reset sector 
366         oT  = 1024; // Reset old strip
367         ooT = t;    // Reset base 
368       }
369       else if (fDebug) {
370         if (ooT == t) 
371           Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]", 
372                d, r, s, t);
373         else 
374           Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
375                "length=%3d (now in %3d, previous %3d)", 
376                d, r, s, ooT, nT, t, oT);
377       }
378         
379       oD = d;
380       oR = r;
381       oS = s;
382       oT = t;
383       nT = TMath::Abs(t-ooT)+1;
384
385       // The longest passage is determined through the angle 
386       Double_t ang  = GetTrackRefTheta(ref, vz);
387       if (ang > angle) {
388         longest = iTrRef;
389         angle   = ang;
390       }
391       // oTheta = ang;
392     } // Loop over track references
393     if (longest < 0) continue; // Nothing found
394
395     // Get the reference corresponding to the longest path through the detector
396     if (fDebug) 
397       Info("Process", "I=%3d L=%3d nT=%3d (out of %3d)", 
398            iTr, longest, nT, fMaxConsequtiveStrips);
399     const AliMCParticle* mother = GetMother(iTr, event);
400
401 #ifdef USE_FLOW_WEIGHTS
402     phi = (mother ? mother->Phi() : particle->Phi());
403     eta = (mother ? mother->Eta() : particle->Eta());
404     Double_t pt = (mother ? mother->Pt() : particle->Pt());
405     Int_t    id = (mother ? mother->PdgCode() : 2212);
406         Double_t weight = CalculateWeight(eta, pt, b, phi, rp, id);
407 #else
408         Double_t weight = 1; // // 
409 #endif
410     StoreParticle(particle, mother, longest, vz, nC, nT, output, weight);
411   } // Loop over tracks
412   return kTRUE;
413 }
414
415 //____________________________________________________________________
416 Double_t
417 AliFMDMCTrackDensity::CalculateWeight(Double_t eta, Double_t pt, Double_t b, 
418                                       Double_t phi, Double_t rp, Int_t id) const
419 {
420   static TF1 gaus = TF1("gaus", "gaus", -6, 6);
421   gaus.SetParameters(0.1, 0., 9);
422   //  gaus.SetParameters(0.1, 0., 3);
423   //  gaus.SetParameters(0.1, 0., 15);
424   
425   const Double_t xCumulant2nd4050ALICE[] = {0.00, 0.25, 0.350,
426                                             0.45, 0.55, 0.650, 
427                                             0.75, 0.85, 0.950,
428                                             1.10, 1.30, 1.500,
429                                             1.70, 1.90, 2.250,
430                                             2.75, 3.25, 3.750,
431                                             4.50};
432   const Double_t yCumulant2nd4050ALICE[] = {0.00000, 0.043400,
433                                             0.059911,0.073516,
434                                             0.089756,0.105486,
435                                             0.117391,0.128199,
436                                             0.138013,0.158271,
437                                             0.177726,0.196383,
438                                             0.208277,0.216648,
439                                             0.242954,0.249961,
440                                             0.240131,0.269006,
441                                             0.207796};
442   const Int_t nPointsCumulant2nd4050ALICE = 
443     sizeof(xCumulant2nd4050ALICE)/sizeof(Double_t);                                      
444   static TGraph alicePointsPt2(nPointsCumulant2nd4050ALICE,xCumulant2nd4050ALICE,yCumulant2nd4050ALICE);
445 #if 0
446   const Double_t xCumulant4th3040ALICE[] = {0.00,0.250,0.35,
447                                             0.45,0.550,0.65,
448                                             0.75,0.850,0.95,
449                                             1.10,1.300,1.50,
450                                             1.70,1.900,2.25,
451                                             2.75,3.250,3.75,
452                                             4.50,5.500,7.00,
453                                             9.000000};
454   const Double_t yCumulant4th3040ALICE[] = {0.000000,0.037071,
455                                             0.048566,0.061083,
456                                             0.070910,0.078831,
457                                             0.091396,0.102026,
458                                             0.109691,0.124449,
459                                             0.139819,0.155561,
460                                             0.165701,0.173678,
461                                             0.191149,0.202015,
462                                             0.204540,0.212560,
463                                             0.195885,0.000000,
464                                             0.000000,0.000000};
465 #endif
466   const Double_t xCumulant4th4050ALICE[] = {0.00,0.25,0.350,
467                                             0.45,0.55,0.650,
468                                             0.75,0.85,0.950,
469                                             1.10,1.30,1.500,
470                                             1.70,1.90,2.250,
471                                             2.75,3.25,3.750,
472                                             4.50};
473   const Double_t yCumulant4th4050ALICE[] = {0.000000,0.038646,
474                                             0.049824,0.066662,
475                                             0.075856,0.081583,
476                                             0.099778,0.104674,
477                                             0.118545,0.131874,
478                                             0.152959,0.155348,
479                                             0.169751,0.179052,
480                                             0.178532,0.198851,
481                                             0.185737,0.239901,
482                                             0.186098};
483   const Int_t nPointsCumulant4th4050ALICE = 
484     sizeof(xCumulant4th4050ALICE)/sizeof(Double_t);   
485   static TGraph alicePointsPt4(nPointsCumulant4th4050ALICE, 
486                                xCumulant4th4050ALICE, 
487                                yCumulant4th4050ALICE);
488
489   const Double_t xCumulant4thTPCrefMultTPConlyAll[] = {1.75,
490                                                        4.225,
491                                                        5.965,
492                                                        7.765,
493                                                        9.215,
494                                                        10.46,
495                                                        11.565,
496                                                        12.575};
497   const Double_t yCumulant4thTPCrefMultTPConlyAll[] = {0.017855,0.032440,
498                                                        0.055818,0.073137,
499                                                        0.083898,0.086690,
500                                                        0.082040,0.077777};
501   const Int_t nPointsCumulant4thTPCrefMultTPConlyAll = 
502     sizeof(xCumulant4thTPCrefMultTPConlyAll)/sizeof(Double_t);
503   TGraph aliceCent(nPointsCumulant4thTPCrefMultTPConlyAll,
504                    xCumulant4thTPCrefMultTPConlyAll,
505                    yCumulant4thTPCrefMultTPConlyAll);
506
507
508   Double_t weight = (20. * gaus.Eval(eta) * (alicePointsPt2.Eval(pt) * 0.5 + 
509                                              alicePointsPt4.Eval(pt) * 0.5) 
510                      * (aliceCent.Eval(b) / aliceCent.Eval(10.46)) 
511                      * 2. * TMath::Cos(2. * (phi - rp)));
512   if      (TMath::Abs(id) == 211)  weight *= 1.3; //pion flow
513   else if (TMath::Abs(id) == 2212) weight *= 1.0;  //proton flow
514   else                             weight *= 0.7;
515   
516   return weight;
517 }
518 //____________________________________________________________________
519 void
520 AliFMDMCTrackDensity::Print(Option_t* /*option*/) const 
521 {
522   char ind[gROOT->GetDirLevel()+1];
523   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
524   ind[gROOT->GetDirLevel()] = '\0';
525   std::cout << ind << ClassName() << ": " << GetName() << '\n'
526             << std::boolalpha 
527             << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
528             << ind << " Max cluster size:       " << fMaxConsequtiveStrips
529             << std::noboolalpha << std::endl;
530   
531 }
532
533 //____________________________________________________________________
534 //
535 // EOF
536 //