4db648a173954d99261832d72b7799f08e2c3925
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDCorrections.cxx
1 #include "AliFMDCorrections.h"
2 #include <AliESDFMD.h>
3 #include <TAxis.h>
4 #include <TList.h>
5 #include <TMath.h>
6 #include "AliForwardCorrectionManager.h"
7 // #include "AliFMDAnaParameters.h"
8 #include "AliLog.h"
9 #include <TH2D.h>
10 #include <TROOT.h>
11 #include <iostream>
12 #include <iomanip>
13
14 ClassImp(AliFMDCorrections)
15 #if 0
16 ; // For Emacs
17 #endif 
18
19 //____________________________________________________________________
20 AliFMDCorrections::AliFMDCorrections()
21   : TNamed(), 
22     fRingHistos(),
23     fDebug(0)
24 {}
25
26 //____________________________________________________________________
27 AliFMDCorrections::AliFMDCorrections(const char* title)
28   : TNamed("fmdCorrections", title), 
29     fRingHistos(), 
30     fDebug(0)
31 {
32   fRingHistos.SetName(GetName());
33   fRingHistos.Add(new RingHistos(1, 'I'));
34   fRingHistos.Add(new RingHistos(2, 'I'));
35   fRingHistos.Add(new RingHistos(2, 'O'));
36   fRingHistos.Add(new RingHistos(3, 'I'));
37   fRingHistos.Add(new RingHistos(3, 'O'));
38 }
39
40 //____________________________________________________________________
41 AliFMDCorrections::AliFMDCorrections(const AliFMDCorrections& o)
42   : TNamed(o), 
43     fRingHistos(), 
44     fDebug(o.fDebug)
45 {
46   TIter    next(&o.fRingHistos);
47   TObject* obj = 0;
48   while ((obj = next())) fRingHistos.Add(obj);
49 }
50
51 //____________________________________________________________________
52 AliFMDCorrections::~AliFMDCorrections()
53 {
54   fRingHistos.Delete();
55 }
56
57 //____________________________________________________________________
58 AliFMDCorrections&
59 AliFMDCorrections::operator=(const AliFMDCorrections& o)
60 {
61   TNamed::operator=(o);
62
63   fDebug   = o.fDebug;
64   fRingHistos.Delete();
65   TIter    next(&o.fRingHistos);
66   TObject* obj = 0;
67   while ((obj = next())) fRingHistos.Add(obj);
68   
69   return *this;
70 }
71
72 //____________________________________________________________________
73 AliFMDCorrections::RingHistos*
74 AliFMDCorrections::GetRingHistos(UShort_t d, Char_t r) const
75 {
76   Int_t idx = -1;
77   switch (d) { 
78   case 1: idx = 0; break;
79   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
80   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
81   }
82   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
83   
84   return static_cast<RingHistos*>(fRingHistos.At(idx));
85 }
86     
87 //____________________________________________________________________
88 Bool_t
89 AliFMDCorrections::Correct(AliForwardUtil::Histos& hists,
90                            UShort_t                vtxbin)
91 {
92   // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
93   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
94
95   UShort_t uvb = vtxbin;
96   for (UShort_t d=1; d<=3; d++) { 
97     UShort_t nr = (d == 1 ? 1 : 2);
98     for (UShort_t q=0; q<nr; q++) { 
99       Char_t      r = (q == 0 ? 'I' : 'O');
100       TH2D*       h = hists.Get(d,r);
101       RingHistos* rh= GetRingHistos(d,r);
102       //TH2F*       bg= pars->GetBackgroundCorrection(d, r, vtxbin);
103       //TH2F*       ef= pars->GetEventSelectionEfficiency("INEL",vtxbin,r);
104       TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb);
105       TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
106       if (!bg) { 
107         AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
108                         d, r, uvb));
109         continue;
110       }
111       if (!ef) { 
112         AliWarning(Form("No event vertex bias correction in vertex bin %d",
113                         uvb));
114         continue;
115       }
116
117       // Divide by primary/total ratio
118       h->Divide(bg);
119       
120       // Divide by the event selection efficiency 
121       h->Divide(ef);
122
123       
124       // if(!pars->SharingEffPresent()) { 
125       //   AliWarning("No sharing efficiencies");
126       //   continue;
127       // }
128       // TH1F* sf = pars->GetSharingEfficiencyTrVtx(d,r,vtxbin); 
129       if (!fcm.GetMergingEfficiency()) { 
130         AliWarning("No merging efficiencies");
131         continue;
132       }
133       TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
134       if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) { 
135         AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
136                         d, r, uvb));
137         continue;
138       }
139
140       
141       for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
142         Float_t c  = sf->GetBinContent(ieta);
143         Float_t ec = sf->GetBinError(ieta);
144         
145         if (c == 0) continue;
146
147         for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
148           Double_t m  = h->GetBinContent(ieta, iphi) / c;
149           Double_t em = h->GetBinError(ieta, iphi);
150           
151           Double_t e  = TMath::Sqrt(em * em + (m * ec) * (m * ec));
152           
153           h->SetBinContent(ieta,iphi,m);
154           h->SetBinError(ieta,iphi,e);
155         }
156       }
157
158       rh->fDensity->Add(h);
159     }
160   }
161   
162   return kTRUE;
163 }
164
165 //____________________________________________________________________
166 void
167 AliFMDCorrections::ScaleHistograms(TList* dir, Int_t nEvents)
168 {
169   if (nEvents <= 0) return;
170   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
171   if (!d) return;
172
173   TIter    next(&fRingHistos);
174   RingHistos* o = 0;
175   while ((o = static_cast<RingHistos*>(next())))
176     o->ScaleHistograms(d, nEvents);
177 }
178 //____________________________________________________________________
179 void
180 AliFMDCorrections::DefineOutput(TList* dir)
181 {
182   TList* d = new TList;
183   d->SetName(GetName());
184   dir->Add(d);
185   TIter    next(&fRingHistos);
186   RingHistos* o = 0;
187   while ((o = static_cast<RingHistos*>(next()))) {
188     o->Output(d);
189   }
190 }
191
192 //____________________________________________________________________
193 void
194 AliFMDCorrections::Print(Option_t* /* option */) const
195 {
196   char ind[gROOT->GetDirLevel()+1];
197   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
198   ind[gROOT->GetDirLevel()] = '\0';
199   std::cout << ind << "AliFMDCorrections: " << GetName() <<  std::endl;
200 }
201
202 //====================================================================
203 AliFMDCorrections::RingHistos::RingHistos()
204   : AliForwardUtil::RingHistos(), 
205     fDensity(0)
206 {}
207
208 //____________________________________________________________________
209 AliFMDCorrections::RingHistos::RingHistos(UShort_t d, Char_t r)
210   : AliForwardUtil::RingHistos(d,r), 
211     fDensity(0)
212 {
213   fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r), 
214                       Form("in FMD%d%c", d, r), 
215                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
216                       0, 2*TMath::Pi());
217   fDensity->SetDirectory(0);
218   fDensity->SetXTitle("#eta");
219   fDensity->SetYTitle("#phi [radians]");
220   fDensity->SetZTitle("Primary N_{ch} density");
221 }
222 //____________________________________________________________________
223 AliFMDCorrections::RingHistos::RingHistos(const RingHistos& o)
224   : AliForwardUtil::RingHistos(o), 
225     fDensity(o.fDensity)
226 {}
227
228 //____________________________________________________________________
229 AliFMDCorrections::RingHistos&
230 AliFMDCorrections::RingHistos::operator=(const RingHistos& o)
231 {
232   AliForwardUtil::RingHistos::operator=(o);
233   
234   if (fDensity) delete fDensity;
235   
236   fDensity = static_cast<TH2D*>(o.fDensity->Clone());
237   
238   return *this;
239 }
240 //____________________________________________________________________
241 AliFMDCorrections::RingHistos::~RingHistos()
242 {
243   if (fDensity) delete fDensity;
244 }
245
246 //____________________________________________________________________
247 void
248 AliFMDCorrections::RingHistos::Output(TList* dir)
249 {
250   TList* d = DefineOutputList(dir);
251   d->Add(fDensity);
252 }
253
254 //____________________________________________________________________
255 void
256 AliFMDCorrections::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
257
258   TList* l = GetOutputList(dir);
259   if (!l) return; 
260
261   TH1* density = GetOutputHist(l,Form("%s_Primary_Density", fName.Data()));
262   if (density) density->Scale(1./nEvents);
263 }
264
265 //____________________________________________________________________
266 //
267 // EOF
268 //
269           
270
271