Next round of upgrades and cleanups of code. The code is now more streamlined and...
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnaParameters.cxx
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
25 //#include <AliCDBManager.h>         // ALICDBMANAGER_H
26 //#include <AliCDBEntry.h>           // ALICDBMANAGER_H
27 //#include "AliFMDRing.h"
28 #include <AliLog.h>
29 #include <Riostream.h>
30 #include <sstream>
31 #include <TSystem.h>
32 #include <TH2D.h>
33 #include <TF1.h>
34 #include <TMath.h>
35 #include "AliESDEvent.h"
36 #include "AliESDVertex.h"
37
38 //====================================================================
39 ClassImp(AliFMDAnaParameters)
40 #if 0
41   ; // This is here to keep Emacs for indenting the next line
42 #endif
43
44 //const char* AliFMDAnaParameters::fgkBackgroundCorrection  = "FMD/Correction/Background";
45 //const char* AliFMDAnaParameters::fgkEnergyDists = "FMD/Correction/EnergyDistribution";
46 const char* AliFMDAnaParameters::fgkBackgroundID = "background";
47 const char* AliFMDAnaParameters::fgkEnergyDistributionID = "energydistributions";
48 const char* AliFMDAnaParameters::fgkEventSelectionEffID  = "eventselectionefficiency";
49 //____________________________________________________________________
50 AliFMDAnaParameters* AliFMDAnaParameters::fgInstance = 0;
51
52 //____________________________________________________________________
53
54 AliFMDAnaParameters* 
55 AliFMDAnaParameters::Instance() 
56 {
57   // Get static instance 
58   if (!fgInstance) fgInstance = new AliFMDAnaParameters;
59   return fgInstance;
60 }
61
62 //____________________________________________________________________
63 AliFMDAnaParameters::AliFMDAnaParameters() :
64   fIsInit(kFALSE),
65   fBackground(0),
66   fEnergyDistribution(0),
67   fEventSelectionEfficiency(0),
68   fCorner1(4.2231, 26.6638),
69   fCorner2(1.8357, 27.9500),
70   fEnergyPath("$ALICE_ROOT/FMD/Correction/EnergyDistribution/energydistributions.root"),
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(kMB1)
76 {
77   
78   
79   //fVerticies.Add(new TVector2(4.2231, 26.6638));
80   // fVerticies.Add(new TVector2(1.8357, 27.9500));
81   // Default constructor 
82 }
83 //____________________________________________________________________
84 void 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();
92   if (what & kEventSelectionEfficiency)   InitEventSelectionEff();
93   
94   fIsInit = kTRUE;
95 }
96 //____________________________________________________________________
97
98 void AliFMDAnaParameters::InitBackground() {
99   
100   //AliCDBEntry*   background = GetEntry(fgkBackgroundCorrection);
101   TFile* fin = TFile::Open(fBackgroundPath.Data());
102   
103   if (!fin) return;
104   
105   fBackground = dynamic_cast<AliFMDAnaCalibBackgroundCorrection*>(fin->Get(fgkBackgroundID));
106   if (!fBackground) AliFatal("Invalid background object from CDB");
107   
108 }
109
110 //____________________________________________________________________
111
112 void AliFMDAnaParameters::InitEnergyDists() {
113   
114   TFile* fin = TFile::Open(fEnergyPath.Data());
115   //AliCDBEntry*   edist = GetEntry(fgkEnergyDists);
116   if (!fin) return;
117   
118   fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(fin->Get(fgkEnergyDistributionID));
119   
120   if (!fEnergyDistribution) AliFatal("Invalid background object from CDB");
121   
122 }
123
124 //____________________________________________________________________
125
126 void 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   
136 }
137 //____________________________________________________________________
138 Float_t AliFMDAnaParameters::GetVtxCutZ() {
139   
140   if(!fIsInit) {
141     AliWarning("Not initialized yet. Call Init() to remedy");
142     return -1;
143   }
144   
145   return fBackground->GetVtxCutZ();
146 }
147
148 //____________________________________________________________________
149 Int_t AliFMDAnaParameters::GetNvtxBins() {
150   
151   if(!fIsInit) {
152     AliWarning("Not initialized yet. Call Init() to remedy");
153     return -1;
154   }
155   
156   return fBackground->GetNvtxBins();
157 }
158 //____________________________________________________________________
159 TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
160   
161   return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
162 }
163 //____________________________________________________________________
164 Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
165   
166   if(!fIsInit) {
167     AliWarning("Not initialized yet. Call Init() to remedy");
168     return 0;
169   }
170   
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);
178   return sigma;
179 }
180
181
182 //____________________________________________________________________
183 Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
184   
185   if(!fIsInit) {
186     AliWarning("Not initialized yet. Call Init() to remedy");
187     return 0;
188   }
189   
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     
197   Float_t mpv           = fitFunc->GetParameter(1);
198   return mpv;
199 }
200 //____________________________________________________________________
201 Float_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;
211   Float_t twoMIPweight    = fitFunc->GetParameter(3);
212   
213   
214   
215   if(twoMIPweight < 1e-05)
216     twoMIPweight = 0;
217   
218   return twoMIPweight;
219 }
220 //____________________________________________________________________
221 Float_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;
231   Float_t threeMIPweight    = fitFunc->GetParameter(4);
232   
233   if(threeMIPweight < 1e-05)
234     threeMIPweight = 0;
235   
236   Float_t twoMIPweight    = fitFunc->GetParameter(3);
237   
238   if(twoMIPweight < 1e-05)
239     threeMIPweight = 0;
240     
241   return threeMIPweight;
242 }
243 //____________________________________________________________________
244 Int_t AliFMDAnaParameters::GetNetaBins() {
245   return GetBackgroundCorrection(1,'I',0)->GetNbinsX();
246   
247 }
248 //____________________________________________________________________
249 Float_t AliFMDAnaParameters::GetEtaMin() {
250
251   return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmin();
252
253 //____________________________________________________________________
254 Float_t AliFMDAnaParameters::GetEtaMax() {
255
256 return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmax();
257
258 }
259 //____________________________________________________________________
260
261 TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det, 
262                                                    Char_t ring, 
263                                                    Int_t vtxbin) {
264   
265   if(!fIsInit) {
266     AliWarning("Not initialized yet. Call Init() to remedy");
267     return 0;
268   }
269   
270   
271   
272   if(vtxbin > fBackground->GetNvtxBins()) {
273     AliWarning(Form("No background object for vertex bin %d", vtxbin));
274     return 0;
275   } 
276   
277   return fBackground->GetBgCorrection(det,ring,vtxbin);
278 }
279 //____________________________________________________________________
280
281 TH1F* 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 }
291 //_____________________________________________________________________
292 Float_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 //_____________________________________________________________________
301 Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) const{
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 //_____________________________________________________________________
313 Float_t AliFMDAnaParameters::GetMinR(Char_t ring) const{
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 //_____________________________________________________________________
326 void 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 //_____________________________________________________________________
339 Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) const
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 //_____________________________________________________________________
369 Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) const
370 {
371   // AliFMDRing fmdring(ring);
372   // fmdring.Init();
373   Float_t   rad       = GetMaxR(ring)-GetMinR(ring);
374   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
375   Float_t   segment   = rad / nStrips;
376   Float_t   r         = GetMinR(ring) + segment*strip;
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;
403 }
404
405 //_____________________________________________________________________
406
407 void AliFMDAnaParameters::GetVertex(AliESDEvent* esd, Double_t* vertexXYZ) 
408 {
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;
421   }
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
434 //____________________________________________________________________
435 Bool_t AliFMDAnaParameters::IsEventTriggered(AliESDEvent *esd) {
436   // check if the event was triggered
437   ULong64_t triggerMask = esd->GetTriggerMask();
438   
439   // definitions from p-p.cfg
440   ULong64_t spdFO = (1 << 14);
441   ULong64_t v0left = (1 << 11);
442   ULong64_t v0right = (1 << 12);
443   
444   switch (fTrigger) {
445   case kMB1: {
446     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
447       return kTRUE;
448     break;
449   }
450   case kMB2: {
451     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
452       return kTRUE;
453     break;
454   }
455   case kSPDFASTOR: {
456     if (triggerMask & spdFO)
457       return kTRUE;
458     break;
459   }
460   }//switch
461
462   return kFALSE;
463 }
464
465 //____________________________________________________________________
466 Float_t 
467 AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)  
468 {
469   //AliFMDRing fmdring(ring);
470   // fmdring.Init();
471   
472   Float_t rad        = GetMaxR(ring)-GetMinR(ring);
473   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
474   Float_t segment    = rad / nStrips;
475   
476   //TVector2* corner1  = fmdring.GetVertex(2);  
477   // TVector2* corner2  = fmdring.GetVertex(3);
478   
479   SetCorners(ring);
480   /*
481   std::cout<<GetMaxR(ring)<<"   "<<fmdring.GetMaxR()<<std::endl;
482   std::cout<<GetMinR(ring)<<"   "<<fmdring.GetMinR()<<std::endl;
483   std::cout<<corner1->X()<<"   "<<fCorner1.X()<<std::endl;
484   std::cout<<corner2->X()<<"   "<<fCorner2.X()<<std::endl;
485   std::cout<<corner1->Y()<<"   "<<fCorner1.Y()<<std::endl;
486   std::cout<<corner2->Y()<<"   "<<fCorner2.Y()<<std::endl;*/
487   Float_t slope      = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
488   Float_t constant   = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
489   Float_t radius     = GetMinR(ring) + strip*segment;
490   
491   Float_t d          = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
492   
493   Float_t arclength  = GetBaseStripLength(ring,strip);
494   if(d>0) {
495     
496     Float_t x        = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
497     Float_t y        = slope*x + constant;
498     Float_t theta    = TMath::ATan2(x,y);
499     
500     if(x < fCorner1.X() && y > fCorner1.Y()) {
501       arclength = radius*theta;                        //One sector since theta is by definition half-hybrid
502       
503     }
504     
505   }
506   
507   return arclength;
508   
509   
510 }
511 //____________________________________________________________________
512 Float_t 
513 AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)  
514 {  
515   // AliFMDRing fmdring(ring);
516   // fmdring.Init();
517   Float_t rad             = GetMaxR(ring)-GetMinR(ring);
518   Float_t nStrips         = (ring == 'I' ? 512 : 256);
519   Float_t nSec            = (ring == 'I' ? 20 : 40);
520   Float_t segment         = rad / nStrips;
521   Float_t basearc         = 2*TMath::Pi() / (0.5*nSec); // One hybrid: 36 degrees inner, 18 outer
522   Float_t radius          = GetMinR(ring) + strip*segment;
523   Float_t basearclength   = 0.5*basearc * radius;                // One sector   
524   
525   return basearclength;
526 }
527 //____________________________________________________________________
528 //
529 // EOF
530 //