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