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