1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
20 //Author: Hans Hjersing Dalsgaard, NBI, hans.dalsgaard@cern.ch
23 #include "AliFMDDebug.h" // ALILOG_H
24 #include "AliFMDAnaParameters.h" // ALIFMDPARAMETERS_H
25 //#include <AliCDBManager.h> // ALICDBMANAGER_H
26 //#include <AliCDBEntry.h> // ALICDBMANAGER_H
27 //#include "AliFMDRing.h"
29 #include <Riostream.h>
36 //====================================================================
37 ClassImp(AliFMDAnaParameters)
39 ; // This is here to keep Emacs for indenting the next line
42 const char* AliFMDAnaParameters::fgkBackgroundCorrection = "FMD/Correction/Background";
43 const char* AliFMDAnaParameters::fgkEnergyDists = "FMD/Correction/EnergyDistribution";
44 const char* AliFMDAnaParameters::fkBackgroundID = "background";
45 const char* AliFMDAnaParameters::fkEnergyDistributionID = "energydistributions";
46 //____________________________________________________________________
47 AliFMDAnaParameters* AliFMDAnaParameters::fgInstance = 0;
49 //____________________________________________________________________
52 AliFMDAnaParameters::Instance()
54 // Get static instance
55 if (!fgInstance) fgInstance = new AliFMDAnaParameters;
59 //____________________________________________________________________
60 AliFMDAnaParameters::AliFMDAnaParameters() :
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")
71 //fVerticies.Add(new TVector2(4.2231, 26.6638));
72 // fVerticies.Add(new TVector2(1.8357, 27.9500));
73 // Default constructor
75 //____________________________________________________________________
76 void AliFMDAnaParameters::Init(Bool_t forceReInit, UInt_t what)
78 // Initialize the parameters manager. We need to get stuff from the
80 if (forceReInit) fIsInit = kFALSE;
82 if (what & kBackgroundCorrection) InitBackground();
83 if (what & kEnergyDistributions) InitEnergyDists();
88 //____________________________________________________________________
90 void AliFMDAnaParameters::InitBackground() {
92 //AliCDBEntry* background = GetEntry(fgkBackgroundCorrection);
93 TFile* fin = TFile::Open(fBackgroundPath.Data());
97 fBackground = dynamic_cast<AliFMDAnaCalibBackgroundCorrection*>(fin->Get(fkBackgroundID));
98 if (!fBackground) AliFatal("Invalid background object from CDB");
101 //____________________________________________________________________
103 void AliFMDAnaParameters::InitEnergyDists() {
105 TFile* fin = TFile::Open(fEnergyPath.Data());
106 //AliCDBEntry* edist = GetEntry(fgkEnergyDists);
109 fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(fin->Get(fkEnergyDistributionID));
111 if (!fEnergyDistribution) AliFatal("Invalid background object from CDB");
114 //____________________________________________________________________
115 Float_t AliFMDAnaParameters::GetVtxCutZ() {
118 AliWarning("Not initialized yet. Call Init() to remedy");
122 return fBackground->GetVtxCutZ();
125 //____________________________________________________________________
126 Int_t AliFMDAnaParameters::GetNvtxBins() {
129 AliWarning("Not initialized yet. Call Init() to remedy");
133 return fBackground->GetNvtxBins();
135 //____________________________________________________________________
136 TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
138 return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
140 //____________________________________________________________________
141 Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
144 AliWarning("Not initialized yet. Call Init() to remedy");
148 TH1F* hEnergyDist = GetEnergyDistribution(det,ring, eta);
149 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
151 AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
154 Float_t sigma = fitFunc->GetParameter(2);
159 //____________________________________________________________________
160 Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
163 AliWarning("Not initialized yet. Call Init() to remedy");
167 TH1F* hEnergyDist = GetEnergyDistribution(det,ring,eta);
168 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
170 AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
174 Float_t MPV = fitFunc->GetParameter(1);
177 //____________________________________________________________________
178 Float_t AliFMDAnaParameters::Get2MIPWeight(Int_t det, Char_t ring, Float_t eta) {
181 AliWarning("Not initialized yet. Call Init() to remedy");
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);
192 if(TwoMIPweight < 1e-05)
197 //____________________________________________________________________
198 Float_t AliFMDAnaParameters::Get3MIPWeight(Int_t det, Char_t ring, Float_t eta) {
201 AliWarning("Not initialized yet. Call Init() to remedy");
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);
210 if(ThreeMIPweight < 1e-05)
213 Float_t TwoMIPweight = fitFunc->GetParameter(3);
215 if(TwoMIPweight < 1e-05)
218 return ThreeMIPweight;
220 //____________________________________________________________________
221 Int_t AliFMDAnaParameters::GetNetaBins() {
222 return GetBackgroundCorrection(1,'I',0)->GetNbinsX();
225 //____________________________________________________________________
226 Float_t AliFMDAnaParameters::GetEtaMin() {
228 return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmin();
230 //____________________________________________________________________
231 Float_t AliFMDAnaParameters::GetEtaMax() {
233 return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmax();
236 //____________________________________________________________________
238 TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det,
243 AliWarning("Not initialized yet. Call Init() to remedy");
249 if(vtxbin > fBackground->GetNvtxBins()) {
250 AliWarning(Form("No background object for vertex bin %d", vtxbin));
254 return fBackground->GetBgCorrection(det,ring,vtxbin);
256 //_____________________________________________________________________
257 Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) {
264 AliWarning("Unknown ring - must be I or O!");
268 //_____________________________________________________________________
269 Float_t AliFMDAnaParameters::GetMinR(Char_t ring) {
276 AliWarning("Unknown ring - must be I or O!");
281 //_____________________________________________________________________
282 void AliFMDAnaParameters::SetCorners(Char_t ring) {
285 fCorner1.Set(4.9895, 15.3560);
286 fCorner2.Set(1.8007, 17.2000);
289 fCorner1.Set(4.2231, 26.6638);
290 fCorner2.Set(1.8357, 27.9500);
294 //_____________________________________________________________________
295 Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec)
297 Int_t nsec = (ring == 'I' ? 20 : 40);
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;
310 Float_t step = 2*TMath::Pi() / nsec;
313 phi = basephi - sec*step;
315 phi = basephi + sec*step;
318 phi = phi +2*TMath::Pi();
319 if(phi > 2*TMath::Pi() )
320 phi = phi - 2*TMath::Pi();
324 //_____________________________________________________________________
325 Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx)
327 // AliFMDRing fmdring(ring);
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;
334 Int_t hybrid = sec / 2;
337 if(!(hybrid%2)) z = 320.266; else z = 319.766;
339 if(det == 2 && ring == 'I' ) {
340 if(!(hybrid%2)) z = 83.666; else z = 83.166;
342 if(det == 2 && ring == 'O' ) {
343 if(!(hybrid%2)) z = 74.966; else z = 75.466;
345 if(det == 3 && ring == 'I' ) {
346 if(!(hybrid%2)) z = -63.066; else z = -62.566;
348 if(det == 3 && ring == 'O' ) {
349 if(!(hybrid%2)) z = -74.966; else z = -75.466;
352 //std::cout<<det<<" "<<ring<<" "<<sec<<" "<<hybrid<<" "<<z<<std::endl;
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));
360 //____________________________________________________________________
361 AliCDBEntry* AliFMDAnaParameters::GetEntry(const char* path, Bool_t fatal) const
363 // Get an entry from the CDB or via preprocessor
364 AliCDBEntry* entry = 0;
365 AliCDBManager* cdb = AliCDBManager::Instance();
366 entry = cdb->Get(path);
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__,
380 //____________________________________________________________________
382 AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)
384 //AliFMDRing fmdring(ring);
387 Float_t rad = GetMaxR(ring)-GetMinR(ring);
388 Float_t Nstrips = (ring == 'I' ? 512 : 256);
389 Float_t segment = rad / Nstrips;
391 //TVector2* corner1 = fmdring.GetVertex(2);
392 // TVector2* corner2 = fmdring.GetVertex(3);
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;
406 Float_t d = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
408 Float_t arclength = GetBaseStripLength(ring,strip);
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);
415 if(x < fCorner1.X() && y > fCorner1.Y()) {
416 arclength = radius*theta; //One sector since theta is by definition half-hybrid
426 //____________________________________________________________________
428 AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)
430 // AliFMDRing fmdring(ring);
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
440 return basearclength;
442 //____________________________________________________________________