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