]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG2/FORWARD/analysis2/AliFMDMCSharingFilter.cxx
Various improvements
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDMCSharingFilter.cxx
1 //
2 // Class to do the sharing correction for MC data.
3 //
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: 
13 //    - None
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 // 
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>
30 #include "AliFMDStripIndex.h"
31 #include <AliLog.h>
32 #include <TROOT.h>
33 #include <iostream>
34 #include <iomanip>
35
36 ClassImp(AliFMDMCSharingFilter)
37 #if 0
38 ; // This is for Emacs
39 #endif 
40
41
42 //____________________________________________________________________
43 AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
44   : AliFMDSharingFilter(title), 
45     fFMD1i(0),
46     fFMD2i(0),
47     fFMD2o(0),
48     fFMD3i(0),
49     fFMD3o(0),
50     fSumEta(0)
51 {
52   // 
53   // Constructor 
54   // 
55   // Parameters:
56   //    title Title of object  - not significant 
57   //
58   fFMD1i = new TH2D("FMD1i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
59   fFMD2i = new TH2D("FMD2i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
60   fFMD2o = new TH2D("FMD2o_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
61   fFMD3i = new TH2D("FMD3i_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
62   fFMD3o = new TH2D("FMD3o_corr", "Merged vs MC", 21, -.5, 20.5, 300, 0, 15);
63   fFMD1i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
64   fFMD1i->SetXTitle("Hits (MC)");
65   fFMD2i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
66   fFMD2i->SetXTitle("Hits (MC)");
67   fFMD2o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
68   fFMD2o->SetXTitle("Hits (MC)");
69   fFMD3i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
70   fFMD3i->SetXTitle("Hits (MC)");
71   fFMD3o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
72   fFMD3o->SetXTitle("Hits (MC)");
73   fFMD1i->SetDirectory(0);
74   fFMD2i->SetDirectory(0);
75   fFMD2o->SetDirectory(0);
76   fFMD3i->SetDirectory(0);
77   fFMD3o->SetDirectory(0);
78   fSumEta = new TH1D("mcSumEta", "MC INEL Truth", 200, -4, 6);
79   fSumEta->SetXTitle("#eta");
80   fSumEta->SetYTitle("dN_{ch}/d#eta");
81   fSumEta->SetDirectory(0);
82   fSumEta->Sumw2();
83   fSumEta->SetMarkerColor(kOrange+2);
84   fSumEta->SetMarkerStyle(22);
85   fSumEta->SetFillColor(0);
86   fSumEta->SetFillStyle(0);
87   
88 }
89
90 //____________________________________________________________________
91 AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o)
92   : AliFMDSharingFilter(o), 
93     fFMD1i(o.fFMD1i),
94     fFMD2i(o.fFMD2i),
95     fFMD2o(o.fFMD2o),
96     fFMD3i(o.fFMD3i),
97     fFMD3o(o.fFMD3o),
98     fSumEta(o.fSumEta)
99 {
100   // 
101   // Copy constructor 
102   // 
103   // Parameters:
104   //    o Object to copy from 
105   //
106 }
107
108 //____________________________________________________________________
109 AliFMDMCSharingFilter::~AliFMDMCSharingFilter()
110 {
111   // 
112   // Destructor
113   //
114   if (fFMD1i)  delete fFMD1i;
115   if (fFMD2i)  delete fFMD2i;
116   if (fFMD2o)  delete fFMD2o;
117   if (fFMD3i)  delete fFMD3i;
118   if (fFMD3o)  delete fFMD3o;
119   if (fSumEta) delete fSumEta;
120 }
121
122 //____________________________________________________________________
123 AliFMDMCSharingFilter&
124 AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o)
125 {
126   // 
127   // Assignment operator 
128   // 
129   // Parameters:
130   //    o Object to assign from 
131   // 
132   // Return:
133   //    Reference to this 
134   //
135   AliFMDSharingFilter::operator=(o);
136   return *this;
137 }
138
139 //____________________________________________________________________
140 void
141 AliFMDMCSharingFilter::StoreParticle(UShort_t   d, 
142                                      Char_t     r,
143                                      UShort_t   s, 
144                                      UShort_t   t, 
145                                      AliESDFMD& output) const
146 {
147   // 
148   // Store a particle hit in FMD<i>dr</i>[<i>s,t</i>] in @a output
149   // 
150   // Parameters:
151   //    d       Detector
152   //    r       Ring
153   //    s       Sector
154   //    t       Strip
155   //    output  Output ESD object
156   //
157   Double_t old = output.Multiplicity(d,r,s,t);
158   if (old == AliESDFMD::kInvalidMult) old = 0;
159   output.SetMultiplicity(d,r,s,t,old+1);
160 }
161
162 //____________________________________________________________________
163 Bool_t
164 AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input, 
165                                 const AliMCEvent& event,
166                                 Double_t          vz,
167                                 AliESDFMD&        output, 
168                                 TH2D*             primary)
169 {
170   // 
171   // Filter the input kinematics and track references, using 
172   // some of the ESD information
173   // 
174   // Parameters:
175   //    input   Input ESD event
176   //    event   Input MC event
177   //    vz      Vertex position 
178   //    output  Output ESD-like object
179   //    primary Per-event histogram of primaries 
180   //
181   // Return:
182   //    True on succes, false otherwise 
183   //
184   output.Clear();
185
186   // Copy eta values to output 
187   for (UShort_t ed = 1; ed <= 3; ed++) { 
188     UShort_t nq = (ed == 1 ? 1 : 2);
189     for (UShort_t eq = 0; eq < nq; eq++) {
190       Char_t   er = (eq == 0 ? 'I' : 'O');
191       UShort_t ns = (eq == 0 ?  20 :  40);
192       UShort_t nt = (eq == 0 ? 512 : 256);
193       for (UShort_t es = 0; es < ns; es++) 
194         for (UShort_t et = 0; et < nt; et++) 
195           output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
196     }
197   }
198   AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
199   Int_t nTracks   = event.GetNumberOfTracks();
200   for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
201     AliMCParticle* particle = 
202       static_cast<AliMCParticle*>(event.GetTrack(iTr));
203     
204     // Check the returned particle 
205     if (!particle) continue;
206     
207     // Check if this charged and a primary 
208     Bool_t isCharged = particle->Charge() != 0;
209     if (!isCharged) continue;
210     Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
211
212     // Fill 'dn/deta' histogram 
213     if (isPrimary) { 
214       fSumEta->Fill(particle->Eta());
215       primary->Fill(particle->Eta(), particle->Phi());
216     }
217
218     Int_t    nTrRef  = particle->GetNumberOfTrackReferences();
219     Int_t    longest = -1;
220     Double_t angle   = 0;
221     UShort_t oD = 0, oS = 1024, oT = 1024;
222     Char_t   oR = '\0';
223     for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
224       AliTrackReference* ref = particle->GetTrackReference(iTrRef);
225       
226       // Check existence 
227       if (!ref) continue;
228
229       // Check that we hit an FMD element 
230       if (ref->DetectorId() != AliTrackReference::kFMD) 
231         continue;
232
233       // Get the detector coordinates 
234       UShort_t d, s, t;
235       Char_t r;
236       AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
237       if (oD > 0 && oR != '\0' && (d != oD || r != oR)) {
238         longest = -1;
239         StoreParticle(oD, oR, oS, oT, output);
240       }
241
242       oD = d;
243       oR = r;
244       oS = s;
245       oT = t;
246
247       // The longest passage is determined through the angle 
248       Double_t x    = ref->X();
249       Double_t y    = ref->Y();
250       Double_t z    = ref->Z()-vz;
251       Double_t rr   = TMath::Sqrt(x*x+y*y);
252       Double_t theta= TMath::ATan2(rr,z);
253       Double_t ang  = TMath::Abs(TMath::Pi()-theta);
254       if (ang > angle) {
255         longest = iTrRef;
256         angle   = ang;
257       }
258     } // Loop over track references
259     if (longest < 0) continue;
260
261     // Get the reference corresponding to the longest path through the detector
262     AliTrackReference* ref = particle->GetTrackReference(longest);
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     StoreParticle(d,r,s,t,output);
270   } // Loop over tracks
271   return kTRUE;
272 }
273
274 //____________________________________________________________________
275 void
276 AliFMDMCSharingFilter::CompareResults(const AliESDFMD&  esd, 
277                                       const AliESDFMD&  mc)
278 {
279   // 
280   // Compare the result of merging to the monte-carlo truth.  This
281   // fills the correlation histograms
282   // 
283   // Parameters:
284   //    esd  ESD after sharing correction
285   //    mc   MC ESD 
286   //
287
288   // Copy eta values to output 
289   for (UShort_t d = 1; d <= 3; d++) { 
290     UShort_t nq = (d == 1 ? 1 : 2);
291     for (UShort_t q = 0; q < nq; q++) {
292       Char_t   r  = (q == 0 ? 'I' : 'O');
293       UShort_t ns = (q == 0 ?  20 :  40);
294       UShort_t nt = (q == 0 ? 512 : 256);
295       TH2*     co = 0;
296       switch (d) { 
297       case 1: co = fFMD1i; break;
298       case 2: co = (q == 0 ? fFMD2i : fFMD2o); break;
299       case 3: co = (q == 0 ? fFMD3i : fFMD3o); break;
300       }
301
302       for (UShort_t s = 0; s < ns; s++) {
303         for (UShort_t t = 0; t < nt; t++) { 
304           Float_t mEsd = esd.Multiplicity(d, r, s, t);
305           Float_t mMc  = mc.Multiplicity(d, r, s, t);
306
307           co->Fill(mMc, mEsd);
308         } 
309       }
310     }
311   }
312 }
313   
314 //____________________________________________________________________
315 void
316 AliFMDMCSharingFilter::DefineOutput(TList* dir)
317 {
318   // 
319   // Define the output histograms.  These are put in a sub list of the
320   // passed list.   The histograms are merged before the parent task calls 
321   // AliAnalysisTaskSE::Terminate 
322   // 
323   // Parameters:
324   //    dir Directory to add to 
325   //
326   AliFMDSharingFilter::DefineOutput(dir);
327   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
328   TList* cd = new TList;
329   cd->SetName("esd_mc_comparion");
330   d->Add(cd);
331   cd->Add(fFMD1i);
332   cd->Add(fFMD2i);
333   cd->Add(fFMD2o);
334   cd->Add(fFMD3i);
335   cd->Add(fFMD3o);
336   dir->Add(fSumEta);
337 }
338
339 //____________________________________________________________________
340 void
341 AliFMDMCSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
342 {
343   // 
344   // Scale the histograms to the total number of events 
345   // 
346   // Parameters:
347   //    dir     Where the output is 
348   //    nEvents Number of events 
349   //
350   AliFMDSharingFilter::ScaleHistograms(dir, nEvents);
351   TH1D* sumEta = static_cast<TH1D*>(dir->FindObject("mcSumEta"));
352   if (!sumEta) { 
353     AliWarning(Form("No mcSumEta histogram found in output list"));
354     return;
355   }
356   sumEta->Scale(1. / nEvents, "width");
357 }
358
359 //____________________________________________________________________
360 void
361 AliFMDMCSharingFilter::Print(Option_t* option) const
362 {
363   // 
364   // Print information
365   // 
366   // Parameters:
367   //    option Not used 
368   //
369   char ind[gROOT->GetDirLevel()+1];
370   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
371   ind[gROOT->GetDirLevel()] = '\0';
372   AliFMDSharingFilter::Print(option);
373 }
374
375 //____________________________________________________________________
376 //
377 // EOF
378 //