779241ce73a65fd98fe42f9f97d2c58b605dc914
[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   fSpecies(kPP)
79 {
80   
81   
82   //fVerticies.Add(new TVector2(4.2231, 26.6638));
83   // fVerticies.Add(new TVector2(1.8357, 27.9500));
84   // Default constructor 
85 }
86 //____________________________________________________________________
87 char* AliFMDAnaParameters::GetPath(const char* species) {
88   
89   char* path ;
90   
91   if(species == fgkBackgroundID)
92     path = Form("%s/%s_%d_%d_%d.root",
93                 fBackgroundPath.Data(),
94                 fgkBackgroundID,
95                 fEnergy,
96                 fTrigger,
97                 fMagField,
98                 fSpecies,
99                 0,
100                 0);
101   if(species == fgkEnergyDistributionID)
102     path = Form("%s/%s_%d_%d_%d.root",
103                 fEnergyPath.Data(),
104                 fgkEnergyDistributionID,
105                 fEnergy,
106                 fTrigger,
107                 fMagField,
108                 fSpecies,
109                 0,
110                 0);
111   if(species == fgkEventSelectionEffID)
112     path = Form("%s/%s_%d_%d_%d.root",
113                 fEventSelectionEffPath.Data(),
114                 fgkEventSelectionEffID,
115                 fEnergy,
116                 fTrigger,
117                 fMagField,
118                 fSpecies,
119                 0,
120                 0);
121
122   return path;
123 }
124 //____________________________________________________________________
125 void 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();
133   if (what & kEventSelectionEfficiency)   InitEventSelectionEff();
134   
135   fIsInit = kTRUE;
136 }
137 //____________________________________________________________________
138
139 void AliFMDAnaParameters::InitBackground() {
140   
141   //AliCDBEntry*   background = GetEntry(fgkBackgroundCorrection);
142   
143   TFile* fin = TFile::Open(GetPath(fgkBackgroundID));
144   
145   if (!fin) return;
146   
147   fBackground = dynamic_cast<AliFMDAnaCalibBackgroundCorrection*>(fin->Get(fgkBackgroundID));
148   if (!fBackground) AliFatal("Invalid background object from CDB");
149   
150 }
151
152 //____________________________________________________________________
153
154 void AliFMDAnaParameters::InitEnergyDists() {
155   
156   TFile* fin = TFile::Open(GetPath(fgkEnergyDistributionID));
157   //AliCDBEntry*   edist = GetEntry(fgkEnergyDists);
158   if (!fin) return;
159   
160   fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(fin->Get(fgkEnergyDistributionID));
161   
162   if (!fEnergyDistribution) AliFatal("Invalid background object from CDB");
163   
164 }
165
166 //____________________________________________________________________
167
168 void AliFMDAnaParameters::InitEventSelectionEff() {
169   
170   //AliCDBEntry*   background = GetEntry(fgkBackgroundCorrection);
171   TFile* fin = TFile::Open(GetPath(fgkEventSelectionEffID));
172                             
173   if (!fin) return;
174   
175   fEventSelectionEfficiency = dynamic_cast<AliFMDAnaCalibEventSelectionEfficiency*>(fin->Get(fgkEventSelectionEffID));
176   if (!fEventSelectionEfficiency) AliFatal("Invalid background object from CDB");
177   
178 }
179 //____________________________________________________________________
180 Float_t AliFMDAnaParameters::GetVtxCutZ() {
181   
182   if(!fIsInit) {
183     AliWarning("Not initialized yet. Call Init() to remedy");
184     return -1;
185   }
186   
187   return fBackground->GetVtxCutZ();
188 }
189
190 //____________________________________________________________________
191 Int_t AliFMDAnaParameters::GetNvtxBins() {
192   
193   if(!fIsInit) {
194     AliWarning("Not initialized yet. Call Init() to remedy");
195     return -1;
196   }
197   
198   return fBackground->GetNvtxBins();
199 }
200 //____________________________________________________________________
201 TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
202   
203   return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
204 }
205 //____________________________________________________________________
206 Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
207   
208   if(!fIsInit) {
209     AliWarning("Not initialized yet. Call Init() to remedy");
210     return 0;
211   }
212   
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);
220   return sigma;
221 }
222
223
224 //____________________________________________________________________
225 Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
226   
227   if(!fIsInit) {
228     AliWarning("Not initialized yet. Call Init() to remedy");
229     return 0;
230   }
231   
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     
239   Float_t mpv           = fitFunc->GetParameter(1);
240   return mpv;
241 }
242 //____________________________________________________________________
243 Float_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;
253   Float_t twoMIPweight    = fitFunc->GetParameter(3);
254   
255   
256   
257   if(twoMIPweight < 1e-05)
258     twoMIPweight = 0;
259   
260   return twoMIPweight;
261 }
262 //____________________________________________________________________
263 Float_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;
273   Float_t threeMIPweight    = fitFunc->GetParameter(4);
274   
275   if(threeMIPweight < 1e-05)
276     threeMIPweight = 0;
277   
278   Float_t twoMIPweight    = fitFunc->GetParameter(3);
279   
280   if(twoMIPweight < 1e-05)
281     threeMIPweight = 0;
282     
283   return threeMIPweight;
284 }
285 //____________________________________________________________________
286 Int_t AliFMDAnaParameters::GetNetaBins() {
287   return GetBackgroundCorrection(1,'I',0)->GetNbinsX();
288   
289 }
290 //____________________________________________________________________
291 Float_t AliFMDAnaParameters::GetEtaMin() {
292
293   return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmin();
294
295 //____________________________________________________________________
296 Float_t AliFMDAnaParameters::GetEtaMax() {
297
298 return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmax();
299
300 }
301 //____________________________________________________________________
302
303 TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det, 
304                                                    Char_t ring, 
305                                                    Int_t vtxbin) {
306   
307   if(!fIsInit) {
308     AliWarning("Not initialized yet. Call Init() to remedy");
309     return 0;
310   }
311   
312   
313   
314   if(vtxbin > fBackground->GetNvtxBins()) {
315     AliWarning(Form("No background object for vertex bin %d", vtxbin));
316     return 0;
317   } 
318   
319   return fBackground->GetBgCorrection(det,ring,vtxbin);
320 }
321 //____________________________________________________________________
322
323 TH1F* 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 }
333 //_____________________________________________________________________
334 Float_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 //_____________________________________________________________________
343 Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) const{
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 //_____________________________________________________________________
355 Float_t AliFMDAnaParameters::GetMinR(Char_t ring) const{
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 //_____________________________________________________________________
368 void 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 //_____________________________________________________________________
381 Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) const
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 //_____________________________________________________________________
411 Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) const
412 {
413   // AliFMDRing fmdring(ring);
414   // fmdring.Init();
415   Float_t   rad       = GetMaxR(ring)-GetMinR(ring);
416   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
417   Float_t   segment   = rad / nStrips;
418   Float_t   r         = GetMinR(ring) + segment*strip;
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;
445 }
446
447 //_____________________________________________________________________
448
449 void AliFMDAnaParameters::GetVertex(AliESDEvent* esd, Double_t* vertexXYZ) 
450 {
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;
463   }
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 }
475
476 //____________________________________________________________________
477 Bool_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
503
504   return kFALSE;
505 }
506
507 //____________________________________________________________________
508 Float_t 
509 AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)  
510 {
511   //AliFMDRing fmdring(ring);
512   // fmdring.Init();
513   
514   Float_t rad        = GetMaxR(ring)-GetMinR(ring);
515   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
516   Float_t segment    = rad / nStrips;
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 //____________________________________________________________________
554 Float_t 
555 AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)  
556 {  
557   // AliFMDRing fmdring(ring);
558   // fmdring.Init();
559   Float_t rad             = GetMaxR(ring)-GetMinR(ring);
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
564   Float_t radius          = GetMinR(ring) + strip*segment;
565   Float_t basearclength   = 0.5*basearc * radius;                // One sector   
566   
567   return basearclength;
568 }
569 //____________________________________________________________________
570 //
571 // EOF
572 //