added lines for cmake
[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}
f58a4769 256//_____________________________________________________________________
78f6f750 257Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) {
258 Float_t radius = 0;
259 if(ring == 'I')
260 radius = 17.2;
261 else if(ring == 'O')
262 radius = 28.0;
263 else
264 AliWarning("Unknown ring - must be I or O!");
265
266 return radius;
267}
268//_____________________________________________________________________
269Float_t AliFMDAnaParameters::GetMinR(Char_t ring) {
270 Float_t radius = 0;
271 if(ring == 'I')
272 radius = 4.5213;
273 else if(ring == 'O')
274 radius = 15.4;
275 else
276 AliWarning("Unknown ring - must be I or O!");
277
278 return radius;
279
280}
281//_____________________________________________________________________
282void AliFMDAnaParameters::SetCorners(Char_t ring) {
283
284 if(ring == 'I') {
285 fCorner1.Set(4.9895, 15.3560);
286 fCorner2.Set(1.8007, 17.2000);
287 }
288 else {
289 fCorner1.Set(4.2231, 26.6638);
290 fCorner2.Set(1.8357, 27.9500);
291 }
292
293}
294//_____________________________________________________________________
f58a4769 295Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec)
296{
297 Int_t nsec = (ring == 'I' ? 20 : 40);
298 Float_t basephi = 0;
299 if(det == 1)
300 basephi = 1.72787594;
301 if(det == 2 && ring == 'I')
302 basephi = 0.15707963;
303 if(det == 2 && ring == 'O')
304 basephi = 0.078539818;
305 if(det == 3 && ring == 'I')
306 basephi = 2.984513044;
307 if(det == 3 && ring == 'O')
308 basephi = 3.06305289;
309
310 Float_t step = 2*TMath::Pi() / nsec;
311 Float_t phi = 0;
312 if(det == 3)
313 phi = basephi - sec*step;
314 else
315 phi = basephi + sec*step;
316
317 if(phi < 0)
318 phi = phi +2*TMath::Pi();
319 if(phi > 2*TMath::Pi() )
320 phi = phi - 2*TMath::Pi();
321
322 return phi;
323}
324//_____________________________________________________________________
325Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx)
326{
78f6f750 327 // AliFMDRing fmdring(ring);
328 // fmdring.Init();
329 Float_t rad = GetMaxR(ring)-GetMinR(ring);
330 Float_t Nstrips = (ring == 'I' ? 512 : 256);
331 Float_t segment = rad / Nstrips;
332 Float_t r = GetMinR(ring) + segment*strip;
f58a4769 333 Float_t z = 0;
334 Int_t hybrid = sec / 2;
335
336 if(det == 1) {
337 if(!(hybrid%2)) z = 320.266; else z = 319.766;
338 }
339 if(det == 2 && ring == 'I' ) {
340 if(!(hybrid%2)) z = 83.666; else z = 83.166;
341 }
342 if(det == 2 && ring == 'O' ) {
343 if(!(hybrid%2)) z = 74.966; else z = 75.466;
344 }
345 if(det == 3 && ring == 'I' ) {
346 if(!(hybrid%2)) z = -63.066; else z = -62.566;
347 }
348 if(det == 3 && ring == 'O' ) {
349 if(!(hybrid%2)) z = -74.966; else z = -75.466;
350 }
351
352 //std::cout<<det<<" "<<ring<<" "<<sec<<" "<<hybrid<<" "<<z<<std::endl;
353
354 // Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
355 Float_t theta = TMath::ATan2(r,z-zvtx);
356 Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
357
358 return eta;
78f6f750 359}/*
d7346eed 360//____________________________________________________________________
361AliCDBEntry* AliFMDAnaParameters::GetEntry(const char* path, Bool_t fatal) const
362{
363 // Get an entry from the CDB or via preprocessor
364 AliCDBEntry* entry = 0;
365 AliCDBManager* cdb = AliCDBManager::Instance();
366 entry = cdb->Get(path);
367
368 if (!entry) {
369 TString msg(Form("No %s found in CDB, perhaps you need to "
370 "use AliFMDCalibFaker?", path));
371 if (fatal) { AliFatal(msg.Data()); }
372 else AliLog::Message(AliLog::kWarning, msg.Data(), "FMD",
373 "AliFMDParameters", "GetEntry", __FILE__,
374 __LINE__);
375 return 0;
376 }
377 return entry;
378}
78f6f750 379 */
380//____________________________________________________________________
381Float_t
382AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)
383{
384 //AliFMDRing fmdring(ring);
385 // fmdring.Init();
386
387 Float_t rad = GetMaxR(ring)-GetMinR(ring);
388 Float_t Nstrips = (ring == 'I' ? 512 : 256);
389 Float_t segment = rad / Nstrips;
390
391 //TVector2* corner1 = fmdring.GetVertex(2);
392 // TVector2* corner2 = fmdring.GetVertex(3);
393
394 SetCorners(ring);
395 /*
396 std::cout<<GetMaxR(ring)<<" "<<fmdring.GetMaxR()<<std::endl;
397 std::cout<<GetMinR(ring)<<" "<<fmdring.GetMinR()<<std::endl;
398 std::cout<<corner1->X()<<" "<<fCorner1.X()<<std::endl;
399 std::cout<<corner2->X()<<" "<<fCorner2.X()<<std::endl;
400 std::cout<<corner1->Y()<<" "<<fCorner1.Y()<<std::endl;
401 std::cout<<corner2->Y()<<" "<<fCorner2.Y()<<std::endl;*/
402 Float_t slope = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
403 Float_t constant = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
404 Float_t radius = GetMinR(ring) + strip*segment;
405
406 Float_t d = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
407
408 Float_t arclength = GetBaseStripLength(ring,strip);
409 if(d>0) {
410
411 Float_t x = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
412 Float_t y = slope*x + constant;
413 Float_t theta = TMath::ATan2(x,y);
414
415 if(x < fCorner1.X() && y > fCorner1.Y()) {
416 arclength = radius*theta; //One sector since theta is by definition half-hybrid
417
418 }
419
420 }
421
422 return arclength;
423
424
425}
426//____________________________________________________________________
427Float_t
428AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)
429{
430 // AliFMDRing fmdring(ring);
431 // fmdring.Init();
432 Float_t rad = GetMaxR(ring)-GetMinR(ring);
433 Float_t Nstrips = (ring == 'I' ? 512 : 256);
434 Float_t Nsec = (ring == 'I' ? 20 : 40);
435 Float_t segment = rad / Nstrips;
436 Float_t basearc = 2*TMath::Pi() / (0.5*Nsec); // One hybrid: 36 degrees inner, 18 outer
437 Float_t radius = GetMinR(ring) + strip*segment;
438 Float_t basearclength = 0.5*basearc * radius; // One sector
439
440 return basearclength;
441}
d7346eed 442//____________________________________________________________________
443//
444// EOF
445//