Fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDMCDensityCalculator.cxx
1 // 
2 // This class calculates the inclusive charged particle density
3 // in each for the 5 FMD rings based on the MC truth.
4 //
5 // Input:
6 //   - AliMCEvent  MC truth event infromation
7 //
8 // Output:
9 //   - None
10 //
11 // Corrections used: 
12 //   - None
13 //
14 //
15 #include "AliFMDMCDensityCalculator.h"
16 #include <TMath.h>
17 #include "AliForwardCorrectionManager.h"
18 #include "AliFMDStripIndex.h"
19 #include "AliMCEvent.h"
20 #include "AliESDFMD.h"
21 #include "AliLog.h"
22 #include <TH2D.h>
23 #include <TProfile2D.h>
24
25 ClassImp(AliFMDMCDensityCalculator)
26 #if 0
27 ; // For Emacs
28 #endif 
29
30 //____________________________________________________________________
31 AliFMDMCDensityCalculator::~AliFMDMCDensityCalculator()
32 {
33   // 
34   // Destructor 
35   //
36   if (fComps)  fComps->Clear();
37   if (fFMD1i)  delete fFMD1i;
38   if (fFMD2i)  delete fFMD2i;
39   if (fFMD2o)  delete fFMD2o;
40   if (fFMD3i)  delete fFMD3i;
41   if (fFMD3o)  delete fFMD3o;
42   if (fFMD1iC) delete fFMD1iC;
43   if (fFMD2iC) delete fFMD2iC;
44   if (fFMD2oC) delete fFMD2oC;
45   if (fFMD3iC) delete fFMD3iC;
46   if (fFMD3oC) delete fFMD3oC;
47 }
48
49 //____________________________________________________________________
50 AliFMDMCDensityCalculator&
51 AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
52 {
53   // 
54   // Assignement operator
55   // 
56   // Parameters:
57   //    o Object to assign from 
58   // 
59   // Return:
60   //    Reference to this object
61   //
62   AliFMDDensityCalculator::operator=(o);
63   return *this;
64 }
65
66
67 //____________________________________________________________________
68 void
69 AliFMDMCDensityCalculator::Init(const TAxis& eAxis)
70 {
71   // 
72   // Initialize this object 
73   // 
74   // Parameters:
75   //    etaAxis Eta axis to use 
76   //
77   AliFMDDensityCalculator::Init(eAxis);
78   fFMD1i  = Make(1,'I',eAxis);
79   fFMD2i  = Make(2,'I',eAxis);
80   fFMD2o  = Make(2,'O',eAxis);
81   fFMD3i  = Make(3,'I',eAxis);
82   fFMD3o  = Make(3,'O',eAxis); 
83   fFMD1iC = Make(1,'I');
84   fFMD2iC = Make(2,'I');
85   fFMD2oC = Make(2,'O');
86   fFMD3iC = Make(3,'I');
87   fFMD3oC = Make(3,'O');
88  
89   fComps->Add(fFMD1i);
90   fComps->Add(fFMD2i);
91   fComps->Add(fFMD2o);
92   fComps->Add(fFMD3i);
93   fComps->Add(fFMD3o);
94   fComps->Add(fFMD1iC);
95   fComps->Add(fFMD2iC);
96   fComps->Add(fFMD2oC);
97   fComps->Add(fFMD3iC);
98   fComps->Add(fFMD3oC);
99 }
100     
101
102 //____________________________________________________________________
103 Bool_t
104 AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD&        fmd,
105                                        AliForwardUtil::Histos& hists)
106 {
107   // 
108   // Calculate the charged particle density from the MC track references. 
109   // 
110   // Parameters:
111   //    fmd    Forward MC event
112   //    hists  Histograms to fill
113   //    vz     Interaction z coordinate @f$ v_z@f$
114   //    vtxBin bin corresponding to @f$ v_z@f$
115   // 
116   // Return:
117   //    true on success
118   //
119   for (UShort_t d=1; d<=3; d++) { 
120     UShort_t nr = (d == 1 ? 1 : 2);
121     for (UShort_t q=0; q<nr; q++) { 
122       Char_t      r = (q == 0 ? 'I' : 'O');
123       UShort_t    ns= (q == 0 ?  20 :  40);
124       UShort_t    nt= (q == 0 ? 512 : 256);
125       TH2D*       h = hists.Get(d,r);
126
127       for (UShort_t s=0; s<ns; s++) { 
128         for (UShort_t t=0; t<nt; t++) {
129           Float_t mult = fmd.Multiplicity(d,r,s,t);
130           
131           if (mult == 0 || mult > 20) continue;
132
133           Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
134           Float_t eta = fmd.Eta(d,r,s,t);
135
136           Float_t corr = Correction(d, r, s, t, 0, eta, false);
137           Float_t sig  = (corr <= 0 ? 0 : mult / corr);
138           h->Fill(eta,phi,sig);
139         }
140       }
141     }
142   }
143   return kTRUE;
144 }
145
146 //____________________________________________________________________
147 TProfile2D*
148 AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r, 
149                                 const TAxis& axis) const
150 {
151   // 
152   // MAke comparison profiles
153   // 
154   // Parameters:
155   //    d     Detector 
156   //    r     Ring 
157   //    axis  Eta axis 
158   // 
159   // Return:
160   //    Newly allocated profile object
161   //
162   TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
163                                    Form("ESD/MC signal for FMD%d%c", d, r),
164                                    axis.GetNbins(), 
165                                    axis.GetXmin(),
166                                    axis.GetXmax(), 
167                                    (r == 'I' || r == 'i') ? 20 : 40,
168                                    0, 2*TMath::Pi());
169   ret->GetXaxis()->SetTitle("#eta");
170   ret->GetYaxis()->SetTitle("#varphi [degrees]");
171   ret->GetZaxis()->SetTitle("#LT incusive density ESD/MC#GT");
172   ret->SetDirectory(0);
173   return ret;
174 }
175 //____________________________________________________________________
176 TH2D*
177 AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
178 {
179   // 
180   // MAke comparison profiles
181   // 
182   // Parameters:
183   //    d     Detector 
184   //    r     Ring 
185   // 
186   // Return:
187   //    Newly allocated profile object
188   //
189   TH2D* ret = new TH2D(Form("FMD%d%c_corr_mc_esd", d, r),
190                        Form("ESD-MC correlation for FMD%d%c", d, r),
191                        200, 0, 20, 200, 0, 20);
192   ret->GetXaxis()->SetTitle("m_{incl} (MC)");
193   ret->GetYaxis()->SetTitle("#Delta/#Delta_{mp} (ESD)");
194   ret->GetZaxis()->SetTitle("Correlation");
195   ret->SetDirectory(0);
196   return ret;
197 }
198 //____________________________________________________________________
199 void
200 AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
201 {
202   // 
203   // Fill comparison profiles
204   // 
205   // Parameters:
206   //    d    Detector 
207   //    r    Ring 
208   //    esd  ESD histogram
209   //    mc   MC histogram
210   //
211   if (!esd || !mc) return;
212   TProfile2D* p = 0;
213   TH2D*       h = 0;
214   switch (d) { 
215   case 1:  
216     p = fFMD1i;                                   
217     h = fFMD1iC;
218     break;
219   case 2:  
220     p = (r == 'I' || r == 'i' ? fFMD2i  : fFMD2o); 
221     h = (r == 'I' || r == 'i' ? fFMD2iC : fFMD2oC); 
222     break;
223   case 3:  
224     p = (r == 'I' || r == 'i' ? fFMD3i  : fFMD3o); 
225     h = (r == 'I' || r == 'i' ? fFMD3iC : fFMD3oC); 
226     break;
227   }
228   if (!p) return;
229
230   for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) { 
231     Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
232     for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) { 
233       Double_t phi  = esd->GetYaxis()->GetBinCenter(iPhi);
234       Double_t mEsd = esd->GetBinContent(iEta,iPhi);
235       Double_t mMc  = mc->GetBinContent(iEta,iPhi);
236       
237       p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
238       h->Fill(mMc,mEsd);
239     }
240   }
241 }
242
243 //____________________________________________________________________
244 Bool_t
245 AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
246                                           AliForwardUtil::Histos& mc)
247 {
248   // 
249   // Compare the result of analysing the ESD for 
250   // the inclusive charged particle density to analysing 
251   // MC truth 
252   // 
253   // Parameters:
254   //    esd 
255   //    mc 
256   // 
257   // Return:
258   //   true 
259   //
260   Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
261   Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
262   Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
263   Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
264   Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
265
266   return kTRUE;
267 }
268
269 //____________________________________________________________________
270 void
271 AliFMDMCDensityCalculator::DefineOutput(TList* dir)
272 {
273   // 
274   // Output diagnostic histograms to directory 
275   // 
276   // Parameters:
277   //    dir List to write in
278   //  
279   AliFMDDensityCalculator::DefineOutput(dir);
280   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
281
282   fComps = new TList;
283   fComps->SetName("esd_mc_comparison");
284   d->Add(fComps);
285
286 }
287 //____________________________________________________________________
288 //
289 // EOF
290 //
291           
292
293