removed the dependence on PWG0 introduced by the implementation of the trigger
[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>
b64db9b1 35#include "AliESDEvent.h"
36#include "AliESDVertex.h"
d7346eed 37
38//====================================================================
39ClassImp(AliFMDAnaParameters)
40#if 0
41 ; // This is here to keep Emacs for indenting the next line
42#endif
43
d05586f1 44//const char* AliFMDAnaParameters::fgkBackgroundCorrection = "FMD/Correction/Background";
45//const char* AliFMDAnaParameters::fgkEnergyDists = "FMD/Correction/EnergyDistribution";
46const char* AliFMDAnaParameters::fgkBackgroundID = "background";
47const char* AliFMDAnaParameters::fgkEnergyDistributionID = "energydistributions";
b64db9b1 48const char* AliFMDAnaParameters::fgkEventSelectionEffID = "eventselectionefficiency";
d7346eed 49//____________________________________________________________________
50AliFMDAnaParameters* AliFMDAnaParameters::fgInstance = 0;
51
52//____________________________________________________________________
53
54AliFMDAnaParameters*
55AliFMDAnaParameters::Instance()
56{
57 // Get static instance
58 if (!fgInstance) fgInstance = new AliFMDAnaParameters;
59 return fgInstance;
60}
61
62//____________________________________________________________________
63AliFMDAnaParameters::AliFMDAnaParameters() :
64 fIsInit(kFALSE),
b82e76e0 65 fBackground(0),
78f6f750 66 fEnergyDistribution(0),
b64db9b1 67 fEventSelectionEfficiency(0),
78f6f750 68 fCorner1(4.2231, 26.6638),
69 fCorner2(1.8357, 27.9500),
70 fEnergyPath("$ALICE_ROOT/FMD/Correction/EnergyDistribution/energydistributions.root"),
b64db9b1 71 fBackgroundPath("$ALICE_ROOT/FMD/Correction/Background/background.root"),
72 fEventSelectionEffPath("$ALICE_ROOT/FMD/Correction/EventSelectionEfficiency/eventselectionefficiency.root"),
73 fProcessPrimary(kFALSE),
74 fProcessHits(kFALSE),
75 fTrigger(AliPWG0Helper::kMB1)
d7346eed 76{
77
78f6f750 78
79 //fVerticies.Add(new TVector2(4.2231, 26.6638));
80 // fVerticies.Add(new TVector2(1.8357, 27.9500));
d7346eed 81 // Default constructor
82}
83//____________________________________________________________________
84void AliFMDAnaParameters::Init(Bool_t forceReInit, UInt_t what)
85{
86 // Initialize the parameters manager. We need to get stuff from the
87 // CDB here.
88 if (forceReInit) fIsInit = kFALSE;
89 if (fIsInit) return;
90 if (what & kBackgroundCorrection) InitBackground();
91 if (what & kEnergyDistributions) InitEnergyDists();
b64db9b1 92 if (what & kEventSelectionEfficiency) InitEventSelectionEff();
78f6f750 93
d7346eed 94 fIsInit = kTRUE;
95}
96//____________________________________________________________________
97
98void AliFMDAnaParameters::InitBackground() {
99
78f6f750 100 //AliCDBEntry* background = GetEntry(fgkBackgroundCorrection);
101 TFile* fin = TFile::Open(fBackgroundPath.Data());
102
103 if (!fin) return;
d7346eed 104
d05586f1 105 fBackground = dynamic_cast<AliFMDAnaCalibBackgroundCorrection*>(fin->Get(fgkBackgroundID));
b82e76e0 106 if (!fBackground) AliFatal("Invalid background object from CDB");
d7346eed 107
108}
b64db9b1 109
d7346eed 110//____________________________________________________________________
111
112void AliFMDAnaParameters::InitEnergyDists() {
113
78f6f750 114 TFile* fin = TFile::Open(fEnergyPath.Data());
115 //AliCDBEntry* edist = GetEntry(fgkEnergyDists);
116 if (!fin) return;
d7346eed 117
d05586f1 118 fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(fin->Get(fgkEnergyDistributionID));
d7346eed 119
b82e76e0 120 if (!fEnergyDistribution) AliFatal("Invalid background object from CDB");
d7346eed 121
b64db9b1 122}
123
124//____________________________________________________________________
125
126void AliFMDAnaParameters::InitEventSelectionEff() {
127
128 //AliCDBEntry* background = GetEntry(fgkBackgroundCorrection);
129 TFile* fin = TFile::Open(fEventSelectionEffPath.Data());
130
131 if (!fin) return;
132
133 fEventSelectionEfficiency = dynamic_cast<AliFMDAnaCalibEventSelectionEfficiency*>(fin->Get(fgkEventSelectionEffID));
134 if (!fEventSelectionEfficiency) AliFatal("Invalid background object from CDB");
135
d7346eed 136}
137//____________________________________________________________________
138Float_t AliFMDAnaParameters::GetVtxCutZ() {
139
140 if(!fIsInit) {
141 AliWarning("Not initialized yet. Call Init() to remedy");
142 return -1;
143 }
144
b82e76e0 145 return fBackground->GetVtxCutZ();
d7346eed 146}
147
148//____________________________________________________________________
149Int_t AliFMDAnaParameters::GetNvtxBins() {
150
151 if(!fIsInit) {
152 AliWarning("Not initialized yet. Call Init() to remedy");
153 return -1;
154 }
155
b82e76e0 156 return fBackground->GetNvtxBins();
d7346eed 157}
158//____________________________________________________________________
5754671c 159TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
d7346eed 160
5754671c 161 return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
8dc823cc 162}
163//____________________________________________________________________
5754671c 164Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
8dc823cc 165
d7346eed 166 if(!fIsInit) {
167 AliWarning("Not initialized yet. Call Init() to remedy");
168 return 0;
169 }
170
5754671c 171 TH1F* hEnergyDist = GetEnergyDistribution(det,ring, eta);
172 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
173 if(!fitFunc) {
174 AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
175 return 1024;
176 }
177 Float_t sigma = fitFunc->GetParameter(2);
8dc823cc 178 return sigma;
179}
180
181
182//____________________________________________________________________
5754671c 183Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
8dc823cc 184
185 if(!fIsInit) {
186 AliWarning("Not initialized yet. Call Init() to remedy");
187 return 0;
188 }
189
5754671c 190 TH1F* hEnergyDist = GetEnergyDistribution(det,ring,eta);
191 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
192 if(!fitFunc) {
193 AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
194 return 1024;
195 }
196
d05586f1 197 Float_t mpv = fitFunc->GetParameter(1);
198 return mpv;
d7346eed 199}
200//____________________________________________________________________
5754671c 201Float_t AliFMDAnaParameters::Get2MIPWeight(Int_t det, Char_t ring, Float_t eta) {
202
203 if(!fIsInit) {
204 AliWarning("Not initialized yet. Call Init() to remedy");
205 return 0;
206 }
207
208 TH1F* hEnergyDist = GetEnergyDistribution(det,ring,eta);
209 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
210 if(!fitFunc) return 0;
d05586f1 211 Float_t twoMIPweight = fitFunc->GetParameter(3);
5754671c 212
213
214
d05586f1 215 if(twoMIPweight < 1e-05)
216 twoMIPweight = 0;
5754671c 217
d05586f1 218 return twoMIPweight;
5754671c 219}
220//____________________________________________________________________
221Float_t AliFMDAnaParameters::Get3MIPWeight(Int_t det, Char_t ring, Float_t eta) {
222
223 if(!fIsInit) {
224 AliWarning("Not initialized yet. Call Init() to remedy");
225 return 0;
226 }
227
228 TH1F* hEnergyDist = GetEnergyDistribution(det,ring,eta);
229 TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc");
230 if(!fitFunc) return 0;
d05586f1 231 Float_t threeMIPweight = fitFunc->GetParameter(4);
5754671c 232
d05586f1 233 if(threeMIPweight < 1e-05)
234 threeMIPweight = 0;
5754671c 235
d05586f1 236 Float_t twoMIPweight = fitFunc->GetParameter(3);
5754671c 237
d05586f1 238 if(twoMIPweight < 1e-05)
239 threeMIPweight = 0;
5754671c 240
d05586f1 241 return threeMIPweight;
5754671c 242}
243//____________________________________________________________________
244Int_t AliFMDAnaParameters::GetNetaBins() {
245 return GetBackgroundCorrection(1,'I',0)->GetNbinsX();
246
247}
248//____________________________________________________________________
249Float_t AliFMDAnaParameters::GetEtaMin() {
8dc823cc 250
5754671c 251 return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmin();
252}
253//____________________________________________________________________
254Float_t AliFMDAnaParameters::GetEtaMax() {
255
256return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmax();
8dc823cc 257
5754671c 258}
259//____________________________________________________________________
8dc823cc 260
d7346eed 261TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det,
262 Char_t ring,
263 Int_t vtxbin) {
b82e76e0 264
d7346eed 265 if(!fIsInit) {
266 AliWarning("Not initialized yet. Call Init() to remedy");
267 return 0;
268 }
b82e76e0 269
270
271
272 if(vtxbin > fBackground->GetNvtxBins()) {
d7346eed 273 AliWarning(Form("No background object for vertex bin %d", vtxbin));
274 return 0;
275 }
276
b82e76e0 277 return fBackground->GetBgCorrection(det,ring,vtxbin);
d7346eed 278}
4fb49e43 279//____________________________________________________________________
280
281TH1F* AliFMDAnaParameters::GetDoubleHitCorrection(Int_t det,
282 Char_t ring) {
283
284 if(!fIsInit) {
285 AliWarning("Not initialized yet. Call Init() to remedy");
286 return 0;
287 }
288
289 return fBackground->GetDoubleHitCorrection(det,ring);
290}
f58a4769 291//_____________________________________________________________________
b64db9b1 292Float_t AliFMDAnaParameters::GetEventSelectionEfficiency(Int_t vtxbin) {
293 if(!fIsInit) {
294 AliWarning("Not initialized yet. Call Init() to remedy");
295 return 0;
296 }
297 return fEventSelectionEfficiency->GetCorrection(vtxbin);
298
299}
300//_____________________________________________________________________
d05586f1 301Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) const{
78f6f750 302 Float_t radius = 0;
303 if(ring == 'I')
304 radius = 17.2;
305 else if(ring == 'O')
306 radius = 28.0;
307 else
308 AliWarning("Unknown ring - must be I or O!");
309
310 return radius;
311}
312//_____________________________________________________________________
d05586f1 313Float_t AliFMDAnaParameters::GetMinR(Char_t ring) const{
78f6f750 314 Float_t radius = 0;
315 if(ring == 'I')
316 radius = 4.5213;
317 else if(ring == 'O')
318 radius = 15.4;
319 else
320 AliWarning("Unknown ring - must be I or O!");
321
322 return radius;
323
324}
325//_____________________________________________________________________
326void AliFMDAnaParameters::SetCorners(Char_t ring) {
327
328 if(ring == 'I') {
329 fCorner1.Set(4.9895, 15.3560);
330 fCorner2.Set(1.8007, 17.2000);
331 }
332 else {
333 fCorner1.Set(4.2231, 26.6638);
334 fCorner2.Set(1.8357, 27.9500);
335 }
336
337}
338//_____________________________________________________________________
d05586f1 339Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) const
f58a4769 340{
341 Int_t nsec = (ring == 'I' ? 20 : 40);
342 Float_t basephi = 0;
343 if(det == 1)
344 basephi = 1.72787594;
345 if(det == 2 && ring == 'I')
346 basephi = 0.15707963;
347 if(det == 2 && ring == 'O')
348 basephi = 0.078539818;
349 if(det == 3 && ring == 'I')
350 basephi = 2.984513044;
351 if(det == 3 && ring == 'O')
352 basephi = 3.06305289;
353
354 Float_t step = 2*TMath::Pi() / nsec;
355 Float_t phi = 0;
356 if(det == 3)
357 phi = basephi - sec*step;
358 else
359 phi = basephi + sec*step;
360
361 if(phi < 0)
362 phi = phi +2*TMath::Pi();
363 if(phi > 2*TMath::Pi() )
364 phi = phi - 2*TMath::Pi();
365
366 return phi;
367}
368//_____________________________________________________________________
d05586f1 369Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) const
f58a4769 370{
78f6f750 371 // AliFMDRing fmdring(ring);
372 // fmdring.Init();
373 Float_t rad = GetMaxR(ring)-GetMinR(ring);
d05586f1 374 Float_t nStrips = (ring == 'I' ? 512 : 256);
375 Float_t segment = rad / nStrips;
78f6f750 376 Float_t r = GetMinR(ring) + segment*strip;
f58a4769 377 Float_t z = 0;
378 Int_t hybrid = sec / 2;
379
380 if(det == 1) {
381 if(!(hybrid%2)) z = 320.266; else z = 319.766;
382 }
383 if(det == 2 && ring == 'I' ) {
384 if(!(hybrid%2)) z = 83.666; else z = 83.166;
385 }
386 if(det == 2 && ring == 'O' ) {
387 if(!(hybrid%2)) z = 74.966; else z = 75.466;
388 }
389 if(det == 3 && ring == 'I' ) {
390 if(!(hybrid%2)) z = -63.066; else z = -62.566;
391 }
392 if(det == 3 && ring == 'O' ) {
393 if(!(hybrid%2)) z = -74.966; else z = -75.466;
394 }
395
396 //std::cout<<det<<" "<<ring<<" "<<sec<<" "<<hybrid<<" "<<z<<std::endl;
397
398 // Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
399 Float_t theta = TMath::ATan2(r,z-zvtx);
400 Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
401
402 return eta;
b64db9b1 403}
404
405//_____________________________________________________________________
406
407void AliFMDAnaParameters::GetVertex(AliESDEvent* esd, Double_t* vertexXYZ)
d7346eed 408{
b64db9b1 409 const AliESDVertex* vertex = 0;
410 vertex = esd->GetPrimaryVertex();
411 if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
412 vertex = esd->GetPrimaryVertexSPD();
413 if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
414 vertex = esd->GetPrimaryVertexTPC();
415 if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
416 vertex = esd->GetVertex();
417 if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
418 vertex->GetXYZ(vertexXYZ);
419 //std::cout<<vertex->GetName()<<" "<< vertex->GetTitle() <<" "<< vertex->GetZv()<<std::endl;
420 return;
d7346eed 421 }
b64db9b1 422 else if (esd->GetESDTZERO()) {
423 vertexXYZ[0] = 0;
424 vertexXYZ[1] = 0;
425 vertexXYZ[2] = esd->GetT0zVertex();
426
427 return;
428 }
429
430 return;
431
432}
433//____________________________________________________________________
434Bool_t AliFMDAnaParameters::IsEventTriggered(AliESDEvent* esd) {
435
436
437 Bool_t trigger = AliPWG0Helper::IsEventTriggered(esd, AliPWG0Helper::kMB1);
438 // Bool_t trigger = AliPWG0Helper::IsEventTriggered(esd, fTrigger);
439 return trigger;
d7346eed 440}
b64db9b1 441
78f6f750 442//____________________________________________________________________
443Float_t
444AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)
445{
446 //AliFMDRing fmdring(ring);
447 // fmdring.Init();
448
449 Float_t rad = GetMaxR(ring)-GetMinR(ring);
d05586f1 450 Float_t nStrips = (ring == 'I' ? 512 : 256);
451 Float_t segment = rad / nStrips;
78f6f750 452
453 //TVector2* corner1 = fmdring.GetVertex(2);
454 // TVector2* corner2 = fmdring.GetVertex(3);
455
456 SetCorners(ring);
457 /*
458 std::cout<<GetMaxR(ring)<<" "<<fmdring.GetMaxR()<<std::endl;
459 std::cout<<GetMinR(ring)<<" "<<fmdring.GetMinR()<<std::endl;
460 std::cout<<corner1->X()<<" "<<fCorner1.X()<<std::endl;
461 std::cout<<corner2->X()<<" "<<fCorner2.X()<<std::endl;
462 std::cout<<corner1->Y()<<" "<<fCorner1.Y()<<std::endl;
463 std::cout<<corner2->Y()<<" "<<fCorner2.Y()<<std::endl;*/
464 Float_t slope = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
465 Float_t constant = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
466 Float_t radius = GetMinR(ring) + strip*segment;
467
468 Float_t d = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
469
470 Float_t arclength = GetBaseStripLength(ring,strip);
471 if(d>0) {
472
473 Float_t x = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
474 Float_t y = slope*x + constant;
475 Float_t theta = TMath::ATan2(x,y);
476
477 if(x < fCorner1.X() && y > fCorner1.Y()) {
478 arclength = radius*theta; //One sector since theta is by definition half-hybrid
479
480 }
481
482 }
483
484 return arclength;
485
486
487}
488//____________________________________________________________________
489Float_t
490AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)
491{
492 // AliFMDRing fmdring(ring);
493 // fmdring.Init();
494 Float_t rad = GetMaxR(ring)-GetMinR(ring);
d05586f1 495 Float_t nStrips = (ring == 'I' ? 512 : 256);
496 Float_t nSec = (ring == 'I' ? 20 : 40);
497 Float_t segment = rad / nStrips;
498 Float_t basearc = 2*TMath::Pi() / (0.5*nSec); // One hybrid: 36 degrees inner, 18 outer
78f6f750 499 Float_t radius = GetMinR(ring) + strip*segment;
500 Float_t basearclength = 0.5*basearc * radius; // One sector
501
502 return basearclength;
503}
d7346eed 504//____________________________________________________________________
505//
506// EOF
507//