Protection against division by zero.
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnaParameters.cxx
CommitLineData
d7346eed 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15//
16//The design of this class is based on the AliFMDParameters class. Its purpose
17//is to hold parameters for the analysis such as background correction and
18//fit functions.
19//
20//Author: Hans Hjersing Dalsgaard, NBI, hans.dalsgaard@cern.ch
21//
22
23#include "AliFMDDebug.h" // ALILOG_H
24#include "AliFMDAnaParameters.h" // ALIFMDPARAMETERS_H
78f6f750 25//#include <AliCDBManager.h> // ALICDBMANAGER_H
26//#include <AliCDBEntry.h> // ALICDBMANAGER_H
27//#include "AliFMDRing.h"
d7346eed 28#include <AliLog.h>
29#include <Riostream.h>
30#include <sstream>
31#include <TSystem.h>
32#include <TH2D.h>
8dc823cc 33#include <TF1.h>
f58a4769 34#include <TMath.h>
d7346eed 35
36//====================================================================
37ClassImp(AliFMDAnaParameters)
38#if 0
39 ; // This is here to keep Emacs for indenting the next line
40#endif
41
b82e76e0 42const char* AliFMDAnaParameters::fgkBackgroundCorrection = "FMD/Correction/Background";
78f6f750 43const char* AliFMDAnaParameters::fgkEnergyDists = "FMD/Correction/EnergyDistribution";
44const char* AliFMDAnaParameters::fkBackgroundID = "background";
45const char* AliFMDAnaParameters::fkEnergyDistributionID = "energydistributions";
d7346eed 46//____________________________________________________________________
47AliFMDAnaParameters* AliFMDAnaParameters::fgInstance = 0;
48
49//____________________________________________________________________
50
51AliFMDAnaParameters*
52AliFMDAnaParameters::Instance()
53{
54 // Get static instance
55 if (!fgInstance) fgInstance = new AliFMDAnaParameters;
56 return fgInstance;
57}
58
59//____________________________________________________________________
60AliFMDAnaParameters::AliFMDAnaParameters() :
61 fIsInit(kFALSE),
b82e76e0 62 fBackground(0),
78f6f750 63 fEnergyDistribution(0),
64 fCorner1(4.2231, 26.6638),
65 fCorner2(1.8357, 27.9500),
66 fEnergyPath("$ALICE_ROOT/FMD/Correction/EnergyDistribution/energydistributions.root"),
67 fBackgroundPath("$ALICE_ROOT/FMD/Correction/Background/background.root")
d7346eed 68{
69
78f6f750 70
71 //fVerticies.Add(new TVector2(4.2231, 26.6638));
72 // fVerticies.Add(new TVector2(1.8357, 27.9500));
d7346eed 73 // Default constructor
74}
75//____________________________________________________________________
76void AliFMDAnaParameters::Init(Bool_t forceReInit, UInt_t what)
77{
78 // Initialize the parameters manager. We need to get stuff from the
79 // CDB here.
80 if (forceReInit) fIsInit = kFALSE;
81 if (fIsInit) return;
82 if (what & kBackgroundCorrection) InitBackground();
83 if (what & kEnergyDistributions) InitEnergyDists();
84
78f6f750 85
d7346eed 86 fIsInit = kTRUE;
87}
88//____________________________________________________________________
89
90void AliFMDAnaParameters::InitBackground() {
91
78f6f750 92 //AliCDBEntry* background = GetEntry(fgkBackgroundCorrection);
93 TFile* fin = TFile::Open(fBackgroundPath.Data());
94
95 if (!fin) return;
d7346eed 96
78f6f750 97 fBackground = dynamic_cast<AliFMDAnaCalibBackgroundCorrection*>(fin->Get(fkBackgroundID));
b82e76e0 98 if (!fBackground) AliFatal("Invalid background object from CDB");
d7346eed 99
100}
101//____________________________________________________________________
102
103void AliFMDAnaParameters::InitEnergyDists() {
104
78f6f750 105 TFile* fin = TFile::Open(fEnergyPath.Data());
106 //AliCDBEntry* edist = GetEntry(fgkEnergyDists);
107 if (!fin) return;
d7346eed 108
78f6f750 109 fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(fin->Get(fkEnergyDistributionID));
d7346eed 110
b82e76e0 111 if (!fEnergyDistribution) AliFatal("Invalid background object from CDB");
d7346eed 112
113}
114//____________________________________________________________________
115Float_t AliFMDAnaParameters::GetVtxCutZ() {
116
117 if(!fIsInit) {
118 AliWarning("Not initialized yet. Call Init() to remedy");
119 return -1;
120 }
121
b82e76e0 122 return fBackground->GetVtxCutZ();
d7346eed 123}
124
125//____________________________________________________________________
126Int_t AliFMDAnaParameters::GetNvtxBins() {
127
128 if(!fIsInit) {
129 AliWarning("Not initialized yet. Call Init() to remedy");
130 return -1;
131 }
132
b82e76e0 133 return fBackground->GetNvtxBins();
d7346eed 134}
135//____________________________________________________________________
5754671c 136TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
d7346eed 137
5754671c 138 return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
8dc823cc 139}
140//____________________________________________________________________
5754671c 141Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
8dc823cc 142
d7346eed 143 if(!fIsInit) {
144 AliWarning("Not initialized yet. Call Init() to remedy");
145 return 0;
146 }
147
5754671c 148 TH1F* hEnergyDist = GetEnergyDistribution(det,ring, eta);
149 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
150 if(!fitFunc) {
151 AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
152 return 1024;
153 }
154 Float_t sigma = fitFunc->GetParameter(2);
8dc823cc 155 return sigma;
156}
157
158
159//____________________________________________________________________
5754671c 160Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
8dc823cc 161
162 if(!fIsInit) {
163 AliWarning("Not initialized yet. Call Init() to remedy");
164 return 0;
165 }
166
5754671c 167 TH1F* hEnergyDist = GetEnergyDistribution(det,ring,eta);
168 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
169 if(!fitFunc) {
170 AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
171 return 1024;
172 }
173
174 Float_t MPV = fitFunc->GetParameter(1);
8dc823cc 175 return MPV;
d7346eed 176}
177//____________________________________________________________________
5754671c 178Float_t AliFMDAnaParameters::Get2MIPWeight(Int_t det, Char_t ring, Float_t eta) {
179
180 if(!fIsInit) {
181 AliWarning("Not initialized yet. Call Init() to remedy");
182 return 0;
183 }
184
185 TH1F* hEnergyDist = GetEnergyDistribution(det,ring,eta);
186 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
187 if(!fitFunc) return 0;
188 Float_t TwoMIPweight = fitFunc->GetParameter(3);
189
190
191
192 if(TwoMIPweight < 1e-05)
193 TwoMIPweight = 0;
194
195 return TwoMIPweight;
196}
197//____________________________________________________________________
198Float_t AliFMDAnaParameters::Get3MIPWeight(Int_t det, Char_t ring, Float_t eta) {
199
200 if(!fIsInit) {
201 AliWarning("Not initialized yet. Call Init() to remedy");
202 return 0;
203 }
204
205 TH1F* hEnergyDist = GetEnergyDistribution(det,ring,eta);
206 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
207 if(!fitFunc) return 0;
208 Float_t ThreeMIPweight = fitFunc->GetParameter(4);
209
210 if(ThreeMIPweight < 1e-05)
211 ThreeMIPweight = 0;
212
213 Float_t TwoMIPweight = fitFunc->GetParameter(3);
214
215 if(TwoMIPweight < 1e-05)
216 ThreeMIPweight = 0;
217
218 return ThreeMIPweight;
219}
220//____________________________________________________________________
221Int_t AliFMDAnaParameters::GetNetaBins() {
222 return GetBackgroundCorrection(1,'I',0)->GetNbinsX();
223
224}
225//____________________________________________________________________
226Float_t AliFMDAnaParameters::GetEtaMin() {
8dc823cc 227
5754671c 228 return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmin();
229}
230//____________________________________________________________________
231Float_t AliFMDAnaParameters::GetEtaMax() {
232
233return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmax();
8dc823cc 234
5754671c 235}
236//____________________________________________________________________
8dc823cc 237
d7346eed 238TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det,
239 Char_t ring,
240 Int_t vtxbin) {
b82e76e0 241
d7346eed 242 if(!fIsInit) {
243 AliWarning("Not initialized yet. Call Init() to remedy");
244 return 0;
245 }
b82e76e0 246
247
248
249 if(vtxbin > fBackground->GetNvtxBins()) {
d7346eed 250 AliWarning(Form("No background object for vertex bin %d", vtxbin));
251 return 0;
252 }
253
b82e76e0 254 return fBackground->GetBgCorrection(det,ring,vtxbin);
d7346eed 255}
4fb49e43 256//____________________________________________________________________
257
258TH1F* AliFMDAnaParameters::GetDoubleHitCorrection(Int_t det,
259 Char_t ring) {
260
261 if(!fIsInit) {
262 AliWarning("Not initialized yet. Call Init() to remedy");
263 return 0;
264 }
265
266 return fBackground->GetDoubleHitCorrection(det,ring);
267}
f58a4769 268//_____________________________________________________________________
78f6f750 269Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) {
270 Float_t radius = 0;
271 if(ring == 'I')
272 radius = 17.2;
273 else if(ring == 'O')
274 radius = 28.0;
275 else
276 AliWarning("Unknown ring - must be I or O!");
277
278 return radius;
279}
280//_____________________________________________________________________
281Float_t AliFMDAnaParameters::GetMinR(Char_t ring) {
282 Float_t radius = 0;
283 if(ring == 'I')
284 radius = 4.5213;
285 else if(ring == 'O')
286 radius = 15.4;
287 else
288 AliWarning("Unknown ring - must be I or O!");
289
290 return radius;
291
292}
293//_____________________________________________________________________
294void AliFMDAnaParameters::SetCorners(Char_t ring) {
295
296 if(ring == 'I') {
297 fCorner1.Set(4.9895, 15.3560);
298 fCorner2.Set(1.8007, 17.2000);
299 }
300 else {
301 fCorner1.Set(4.2231, 26.6638);
302 fCorner2.Set(1.8357, 27.9500);
303 }
304
305}
306//_____________________________________________________________________
f58a4769 307Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec)
308{
309 Int_t nsec = (ring == 'I' ? 20 : 40);
310 Float_t basephi = 0;
311 if(det == 1)
312 basephi = 1.72787594;
313 if(det == 2 && ring == 'I')
314 basephi = 0.15707963;
315 if(det == 2 && ring == 'O')
316 basephi = 0.078539818;
317 if(det == 3 && ring == 'I')
318 basephi = 2.984513044;
319 if(det == 3 && ring == 'O')
320 basephi = 3.06305289;
321
322 Float_t step = 2*TMath::Pi() / nsec;
323 Float_t phi = 0;
324 if(det == 3)
325 phi = basephi - sec*step;
326 else
327 phi = basephi + sec*step;
328
329 if(phi < 0)
330 phi = phi +2*TMath::Pi();
331 if(phi > 2*TMath::Pi() )
332 phi = phi - 2*TMath::Pi();
333
334 return phi;
335}
336//_____________________________________________________________________
337Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx)
338{
78f6f750 339 // AliFMDRing fmdring(ring);
340 // fmdring.Init();
341 Float_t rad = GetMaxR(ring)-GetMinR(ring);
342 Float_t Nstrips = (ring == 'I' ? 512 : 256);
343 Float_t segment = rad / Nstrips;
344 Float_t r = GetMinR(ring) + segment*strip;
f58a4769 345 Float_t z = 0;
346 Int_t hybrid = sec / 2;
347
348 if(det == 1) {
349 if(!(hybrid%2)) z = 320.266; else z = 319.766;
350 }
351 if(det == 2 && ring == 'I' ) {
352 if(!(hybrid%2)) z = 83.666; else z = 83.166;
353 }
354 if(det == 2 && ring == 'O' ) {
355 if(!(hybrid%2)) z = 74.966; else z = 75.466;
356 }
357 if(det == 3 && ring == 'I' ) {
358 if(!(hybrid%2)) z = -63.066; else z = -62.566;
359 }
360 if(det == 3 && ring == 'O' ) {
361 if(!(hybrid%2)) z = -74.966; else z = -75.466;
362 }
363
364 //std::cout<<det<<" "<<ring<<" "<<sec<<" "<<hybrid<<" "<<z<<std::endl;
365
366 // Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
367 Float_t theta = TMath::ATan2(r,z-zvtx);
368 Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
369
370 return eta;
78f6f750 371}/*
d7346eed 372//____________________________________________________________________
373AliCDBEntry* AliFMDAnaParameters::GetEntry(const char* path, Bool_t fatal) const
374{
375 // Get an entry from the CDB or via preprocessor
376 AliCDBEntry* entry = 0;
377 AliCDBManager* cdb = AliCDBManager::Instance();
378 entry = cdb->Get(path);
379
380 if (!entry) {
381 TString msg(Form("No %s found in CDB, perhaps you need to "
382 "use AliFMDCalibFaker?", path));
383 if (fatal) { AliFatal(msg.Data()); }
384 else AliLog::Message(AliLog::kWarning, msg.Data(), "FMD",
385 "AliFMDParameters", "GetEntry", __FILE__,
386 __LINE__);
387 return 0;
388 }
389 return entry;
390}
78f6f750 391 */
392//____________________________________________________________________
393Float_t
394AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)
395{
396 //AliFMDRing fmdring(ring);
397 // fmdring.Init();
398
399 Float_t rad = GetMaxR(ring)-GetMinR(ring);
400 Float_t Nstrips = (ring == 'I' ? 512 : 256);
401 Float_t segment = rad / Nstrips;
402
403 //TVector2* corner1 = fmdring.GetVertex(2);
404 // TVector2* corner2 = fmdring.GetVertex(3);
405
406 SetCorners(ring);
407 /*
408 std::cout<<GetMaxR(ring)<<" "<<fmdring.GetMaxR()<<std::endl;
409 std::cout<<GetMinR(ring)<<" "<<fmdring.GetMinR()<<std::endl;
410 std::cout<<corner1->X()<<" "<<fCorner1.X()<<std::endl;
411 std::cout<<corner2->X()<<" "<<fCorner2.X()<<std::endl;
412 std::cout<<corner1->Y()<<" "<<fCorner1.Y()<<std::endl;
413 std::cout<<corner2->Y()<<" "<<fCorner2.Y()<<std::endl;*/
414 Float_t slope = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
415 Float_t constant = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
416 Float_t radius = GetMinR(ring) + strip*segment;
417
418 Float_t d = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
419
420 Float_t arclength = GetBaseStripLength(ring,strip);
421 if(d>0) {
422
423 Float_t x = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
424 Float_t y = slope*x + constant;
425 Float_t theta = TMath::ATan2(x,y);
426
427 if(x < fCorner1.X() && y > fCorner1.Y()) {
428 arclength = radius*theta; //One sector since theta is by definition half-hybrid
429
430 }
431
432 }
433
434 return arclength;
435
436
437}
438//____________________________________________________________________
439Float_t
440AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)
441{
442 // AliFMDRing fmdring(ring);
443 // fmdring.Init();
444 Float_t rad = GetMaxR(ring)-GetMinR(ring);
445 Float_t Nstrips = (ring == 'I' ? 512 : 256);
446 Float_t Nsec = (ring == 'I' ? 20 : 40);
447 Float_t segment = rad / Nstrips;
448 Float_t basearc = 2*TMath::Pi() / (0.5*Nsec); // One hybrid: 36 degrees inner, 18 outer
449 Float_t radius = GetMinR(ring) + strip*segment;
450 Float_t basearclength = 0.5*basearc * radius; // One sector
451
452 return basearclength;
453}
d7346eed 454//____________________________________________________________________
455//
456// EOF
457//