Various upgrades and bug fixes etc
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / 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 "AliTriggerAnalysis.h"
36 #include "AliPhysicsSelection.h"
37 #include "AliBackgroundSelection.h"
38 #include "AliMultiplicity.h"
39 #include "AliESDEvent.h"
40 #include "AliESDVertex.h"
41
42
43 //====================================================================
44 ClassImp(AliFMDAnaParameters)
45 #if 0
46   ; // This is here to keep Emacs for indenting the next line
47 #endif
48
49 //const char* AliFMDAnaParameters::fgkBackgroundCorrection  = "FMD/Correction/Background";
50 //const char* AliFMDAnaParameters::fgkEnergyDists = "FMD/Correction/EnergyDistribution";
51 const char* AliFMDAnaParameters::fgkBackgroundID         = "background";
52 const char* AliFMDAnaParameters::fgkEnergyDistributionID = "energydistributions";
53 const char* AliFMDAnaParameters::fgkEventSelectionEffID  = "eventselectionefficiency";
54 const char* AliFMDAnaParameters::fgkSharingEffID         = "sharingefficiency";
55 //____________________________________________________________________
56 AliFMDAnaParameters* AliFMDAnaParameters::fgInstance = 0;
57
58 //____________________________________________________________________
59
60 AliFMDAnaParameters* 
61 AliFMDAnaParameters::Instance() 
62 {
63   // Get static instance 
64   if (!fgInstance) fgInstance = new AliFMDAnaParameters;
65   return fgInstance;
66 }
67
68 //____________________________________________________________________
69 AliFMDAnaParameters::AliFMDAnaParameters() :
70   fIsInit(kFALSE),
71   fBackground(0),
72   fEnergyDistribution(0),
73   fEventSelectionEfficiency(0),
74   fSharingEfficiency(0),
75   fCorner1(4.2231, 26.6638),
76   fCorner2(1.8357, 27.9500),
77   fEnergyPath("$ALICE_ROOT/PWG2/FORWARD/corrections/EnergyDistribution"),
78   fBackgroundPath("$ALICE_ROOT/PWG2/FORWARD/corrections/Background"),
79   fEventSelectionEffPath("$ALICE_ROOT/PWG2/FORWARD/corrections/EventSelectionEfficiency"),
80   fSharingEffPath("$ALICE_ROOT/PWG2/FORWARD/corrections/SharingEfficiency"),
81   fProcessPrimary(kFALSE),
82   fProcessHits(kFALSE),
83   fTrigger(kMB1),
84   fEnergy(k900),
85   fMagField(k5G),
86   fSpecies(kPP),
87   fPhysicsSelection(0),
88   fRealData(kFALSE),
89   fSPDlowLimit(0),
90   fSPDhighLimit(999999999),
91   fCentralSelection(kFALSE),
92   fSharingObjectPresent(kTRUE),
93   fNumberOfEtaBinsToCut(1),
94   fEtaLowBinLimits(),
95   fEtaHighBinLimits(),
96   fTriggerInel(kFALSE),
97   fTriggerNSD(kFALSE),
98   fTriggerEmpty(kFALSE),
99   fUseBuiltInNSD(kFALSE),
100   fInelGtZero(kFALSE)
101 {
102   // Default constructor 
103   //  fPhysicsSelection = new AliPhysicsSelection;
104   //  fPhysicsSelection->SetAnalyzeMC(kTRUE); //For the background correction. This is reset in Init if relevant
105   // fPhysicsSelection->SetUseBXNumbers(kFALSE);
106   
107   
108   //  AliBackgroundSelection* backgroundSelection = new AliBackgroundSelection("bg","bg");
109   //  backgroundSelection->Init();
110   //  fPhysicsSelection->AddBackgroundIdentification(backgroundSelection);
111   //fPhysicsSelection->Initialize(104792);
112   // Do not use this - it is only for IO 
113   fgInstance = this;
114   
115 }
116 //____________________________________________________________________
117 char* AliFMDAnaParameters::GetPath(const char* species) {
118   //Get path of object
119   char* path = "";
120   
121   if(species == fgkBackgroundID)
122     path = Form("%s/%s_%d_%d_%d_%d_%d_%d.root",
123                 fBackgroundPath.Data(),
124                 fgkBackgroundID,
125                 fEnergy,
126                 fTrigger,
127                 fMagField,
128                 fSpecies,
129                 fInelGtZero,
130                 0);
131   if(species == fgkEnergyDistributionID)
132     path = Form("%s/%s_%d_%d_%d_%d_%d_%d.root",
133                 fEnergyPath.Data(),
134                 fgkEnergyDistributionID,
135                 fEnergy,
136                 fTrigger,
137                 fMagField,
138                 fSpecies,
139                 fRealData,
140                 0);
141   if(species == fgkEventSelectionEffID)
142     path = Form("%s/%s_%d_%d_%d_%d_%d_%d.root",
143                 fEventSelectionEffPath.Data(),
144                 fgkEventSelectionEffID,
145                 fEnergy,
146                 fTrigger,
147                 fMagField,
148                 fSpecies,
149                 fInelGtZero,
150                 0);
151   if(species == fgkSharingEffID)
152     path = Form("%s/%s_%d_%d_%d_%d_%d_%d.root",
153                 fSharingEffPath.Data(),
154                 fgkSharingEffID,
155                 fEnergy,
156                 fTrigger,
157                 fMagField,
158                 fSpecies,
159                 0,
160                 0);
161
162   return path;
163 }
164 //____________________________________________________________________
165 void AliFMDAnaParameters::Init(Bool_t forceReInit, UInt_t what)
166 {
167   // Initialize the parameters manager.  We need to get stuff from files here.
168   
169   /*  AliPhysicsSelection* test = (AliPhysicsSelection*)((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetEventSelection();
170   
171     if(fPhysicsSelection) {
172     if(!fRealData) {
173       fPhysicsSelection->SetAnalyzeMC(kTRUE);
174     }
175     else fPhysicsSelection->SetAnalyzeMC(kFALSE);
176   }
177   */
178   if (forceReInit) fIsInit = kFALSE;
179   if (fIsInit) return;
180   if (what & kBackgroundCorrection)       InitBackground();
181   if (what & kEnergyDistributions)        InitEnergyDists();
182   if (what & kEventSelectionEfficiency)   InitEventSelectionEff();
183   if (what & kSharingEfficiency)          InitSharingEff();
184   
185   fIsInit = kTRUE;
186   
187   if(fBackground)
188     FindEtaLimits();
189   
190  
191     
192 }
193 //____________________________________________________________________
194
195 void AliFMDAnaParameters::InitBackground() {
196   //Init background correction objects.
197   
198   TFile* fin = TFile::Open(GetPath(fgkBackgroundID));
199   
200   if (!fin) return;
201   
202   fBackground = dynamic_cast<AliFMDAnaCalibBackgroundCorrection*>(fin->Get(fgkBackgroundID));
203   if (!fBackground) AliFatal("Invalid background object from CDB");
204   
205   
206
207 }
208
209 //____________________________________________________________________
210
211 void AliFMDAnaParameters::InitEnergyDists() {
212   //Init energy distributions
213   TFile* fin = TFile::Open(GetPath(fgkEnergyDistributionID));
214   
215   if (!fin) return;
216   
217   fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(fin->Get(fgkEnergyDistributionID));
218   
219   if (!fEnergyDistribution) AliFatal("Invalid background object from CDB");
220   
221 }
222
223 //____________________________________________________________________
224
225 void AliFMDAnaParameters::InitEventSelectionEff() {
226   //Init event selection objects
227
228   TFile* fin = TFile::Open(GetPath(fgkEventSelectionEffID));
229                             
230   if (!fin) return;
231   
232   fEventSelectionEfficiency = dynamic_cast<AliFMDAnaCalibEventSelectionEfficiency*>(fin->Get(fgkEventSelectionEffID));
233   if (!fEventSelectionEfficiency) AliFatal("Invalid background object from CDB");
234   
235 }
236
237 //____________________________________________________________________
238
239 void AliFMDAnaParameters::InitSharingEff() {
240   
241   fSharingObjectPresent = kTRUE;
242   TFile* fin = TFile::Open(GetPath(fgkSharingEffID));
243                             
244   if (!fin) {
245     fSharingObjectPresent = kFALSE;
246     return; 
247   }
248   
249   fSharingEfficiency = dynamic_cast<AliFMDAnaCalibSharingEfficiency*>(fin->Get(fgkSharingEffID));
250   if (!fSharingEfficiency) {
251     fSharingObjectPresent = kFALSE;
252     return; 
253   }
254   
255 }
256 //____________________________________________________________________
257 void AliFMDAnaParameters::FindEtaLimits() {
258
259   fEtaLowBinLimits.SetBins(4,0,4,2,0,2,GetNvtxBins(),0,GetNvtxBins());
260   fEtaHighBinLimits.SetBins(4,0,4,2,0,2,GetNvtxBins(),0,GetNvtxBins());
261   for(Int_t det=0; det<=3;det++) {
262     Int_t nRings = (det<=1 ? 1 : 2);
263     for (UShort_t ir = 0; ir < nRings; ir++) {
264       Char_t ringChar = (ir == 0 ? 'I' : 'O');
265       for(Int_t v =0; v<GetNvtxBins(); v++) {
266         fEtaLowBinLimits.SetBinContent(det,ir,v,GetFirstEtaBinFromMap(v, det, ringChar));
267         fEtaHighBinLimits.SetBinContent(det,ir,v,GetLastEtaBinFromMap(v, det, ringChar));
268         //std::cout<<det<<"   "<<ringChar<<"   "<<fEtaLowBinLimits.GetBinContent(det,ir,v)<<"   "<<fEtaHighBinLimits.GetBinContent(det,ir,v)<<std::endl;
269       }
270     }
271   }
272
273 }
274 //____________________________________________________________________
275
276 void AliFMDAnaParameters::PrintStatus() const
277 {
278   //Print current status
279   TString energystring;
280   switch(fEnergy) {
281   case k900:
282     energystring.Form("900 GeV");   break;
283   case k7000:
284     energystring.Form("7000 GeV");  break;
285   case k10000:
286     energystring.Form("10000 GeV"); break;
287   case k14000:
288     energystring.Form("14000 GeV"); break;
289   default:
290     energystring.Form("invalid energy"); break;
291   }
292   TString triggerstring;
293   switch(fTrigger) {
294   case kMB1:
295     triggerstring.Form("Minimum bias 1");   break;
296   case kMB2:
297     triggerstring.Form("Minimum bias 2");   break;
298   case kSPDFASTOR:
299     triggerstring.Form("SPD FAST OR");   break;
300   case kNOCTP:
301     triggerstring.Form("NO TRIGGER TEST");   break;
302   default:
303     energystring.Form("invalid trigger"); break;
304   }
305   TString magstring;
306   switch(fMagField) {
307   case k5G:
308     magstring.Form("5 kGaus");   break;
309   case k0G:
310     magstring.Form("0 kGaus");   break;
311   default:
312     magstring.Form("invalid mag field"); break;
313   }
314   TString collsystemstring;
315   switch(fSpecies) {
316   case kPP:
317     collsystemstring.Form("p+p");   break;
318   case kPbPb:
319     collsystemstring.Form("Pb+Pb");   break;
320   default:
321     collsystemstring.Form("invalid collision system");   break;
322   }
323   
324   TString datastring;
325   switch(fRealData) {
326   case kTRUE:
327     datastring.Form("Nature"); break;
328   case kFALSE: 
329     datastring.Form("MC"); break;
330   default:
331     datastring.Form("Unknown"); break ;
332   }
333   
334   TString InelString;
335   if(fInelGtZero) InelString = "INEL > 0";
336   else InelString = "INEL";
337   
338   std::cout<<"Energy       = "<<energystring.Data()<<std::endl;
339   std::cout<<"Trigger      = "<<triggerstring.Data()<<std::endl;
340   std::cout<<"Mag Field    = "<<magstring.Data()<<std::endl;
341   std::cout<<"Coll System  = "<<collsystemstring.Data()<<std::endl;
342   std::cout<<"Data origin  = "<<datastring.Data()<<std::endl;
343   std::cout<<"Basic trigger: "<<InelString.Data()<<std::endl;
344   
345   
346 }
347
348 //____________________________________________________________________
349 Float_t AliFMDAnaParameters::GetVtxCutZ() {
350   //Get the z vtx cut in analysis
351   if(!fIsInit) {
352     AliWarning("Not initialized yet. Call Init() to remedy");
353     return -1;
354   }
355   
356   return fBackground->GetVtxCutZ();
357 }
358
359 //____________________________________________________________________
360 Int_t AliFMDAnaParameters::GetNvtxBins() {
361   //Get number of vtx bins
362   if(!fIsInit) {
363     AliWarning("Not initialized yet. Call Init() to remedy");
364     return -1;
365   }
366   
367   return fBackground->GetNvtxBins();
368 }
369 //____________________________________________________________________
370 TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
371   
372   return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
373 }
374 //____________________________________________________________________
375 TH1F* AliFMDAnaParameters::GetEmptyEnergyDistribution(Int_t det, Char_t ring) {
376   
377   return fEnergyDistribution->GetEmptyEnergyDistribution(det, ring);
378 }
379 //____________________________________________________________________
380 TH1F* AliFMDAnaParameters::GetRingEnergyDistribution(Int_t det, Char_t ring) {
381   
382   return fEnergyDistribution->GetRingEnergyDistribution(det, ring);
383 }
384 //____________________________________________________________________
385 Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
386   //Get sigma of Landau fits
387   if(!fIsInit) {
388     AliWarning("Not initialized yet. Call Init() to remedy");
389     return 0;
390   }
391   
392   TH1F* hEnergyDist       = GetEnergyDistribution(det,ring, eta);
393   TF1*  fitFunc           = hEnergyDist->GetFunction("FMDfitFunc");
394   if(!fitFunc) {
395     //AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
396     return 1024;
397   }
398   Float_t sigma           = fitFunc->GetParameter(2);
399   return sigma;
400 }
401
402
403 //____________________________________________________________________
404 Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
405   //Get MPV of landau fits
406   if(!fIsInit) {
407     AliWarning("Not initialized yet. Call Init() to remedy");
408     return 0;
409   }
410   //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
411   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
412   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
413   
414   if(!fitFunc) {
415     AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
416     std::cout<<hEnergyDist->GetEntries()<<"    "<<GetEtaBin(eta)<<std::endl;
417     return 1024;
418   }
419   
420   Float_t mpv           = fitFunc->GetParameter(1);
421   return mpv;
422 }
423 //____________________________________________________________________
424 Float_t AliFMDAnaParameters::GetConstant(Int_t det, Char_t ring, Float_t eta) {
425   //Get constant parameter of Landau fits
426   if(!fIsInit) {
427     AliWarning("Not initialized yet. Call Init() to remedy");
428     return 0;
429   }
430   
431   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
432   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
433   if(!fitFunc) {
434     AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
435     return 0;
436   }
437     
438   Float_t mpv           = fitFunc->GetParameter(0);
439   return mpv;
440 }
441 //____________________________________________________________________
442 Float_t AliFMDAnaParameters::Get2MIPWeight(Int_t det, Char_t ring, Float_t eta) {
443   //Get 2 MIP weights of convoluted Landau fits
444   if(!fIsInit) {
445     AliWarning("Not initialized yet. Call Init() to remedy");
446     return 0;
447   }
448   
449   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
450   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
451   if(!fitFunc) return 0;
452   Float_t twoMIPweight    = fitFunc->GetParameter(3);
453   
454   
455   
456   if(twoMIPweight < 1e-05)
457     twoMIPweight = 0;
458   
459   return twoMIPweight;
460 }
461 //____________________________________________________________________
462 Float_t AliFMDAnaParameters::Get3MIPWeight(Int_t det, Char_t ring, Float_t eta) {
463   //Get 3 MIP weights of convoluted Landau fits
464   if(!fIsInit) {
465     AliWarning("Not initialized yet. Call Init() to remedy");
466     return 0;
467   }
468   
469   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
470   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
471   if(!fitFunc) return 0;
472   Float_t threeMIPweight    = fitFunc->GetParameter(4);
473   
474   if(threeMIPweight < 1e-05)
475     threeMIPweight = 0;
476   
477   Float_t twoMIPweight    = fitFunc->GetParameter(3);
478   
479   if(twoMIPweight < 1e-05)
480     threeMIPweight = 0;
481     
482   return threeMIPweight;
483 }
484 //____________________________________________________________________
485 Int_t AliFMDAnaParameters::GetNetaBins() {
486   return GetBackgroundCorrection(1,'I',5)->GetNbinsX();
487   
488 }
489 //____________________________________________________________________
490 Float_t AliFMDAnaParameters::GetEtaMin() {
491
492   return GetBackgroundCorrection(1,'I',5)->GetXaxis()->GetXmin();
493
494 //____________________________________________________________________
495 Float_t AliFMDAnaParameters::GetEtaMax() {
496
497 return GetBackgroundCorrection(1,'I',5)->GetXaxis()->GetXmax();
498
499 }
500 //____________________________________________________________________
501 Int_t AliFMDAnaParameters::GetEtaBin(Float_t eta) {
502   
503   TAxis testaxis(GetNetaBins(),GetEtaMin(),GetEtaMax());
504   Int_t binnumber = testaxis.FindBin(eta) ;
505   
506   return binnumber;
507
508 }
509 //____________________________________________________________________
510
511 TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det, 
512                                                    Char_t ring, 
513                                                    Int_t vtxbin) {
514   //Get background correction histogram
515   if(!fIsInit) {
516     AliWarning("Not initialized yet. Call Init() to remedy");
517     return 0;
518   }
519   
520   
521   
522   if(vtxbin > fBackground->GetNvtxBins()) {
523     AliWarning(Form("No background object for vertex bin %d", vtxbin));
524     return 0;
525   } 
526   
527   return fBackground->GetBgCorrection(det,ring,vtxbin);
528 }
529 //____________________________________________________________________
530 TH2F* AliFMDAnaParameters::GetBackgroundCorrectionNSD(Int_t det, 
531                                                       Char_t ring, 
532                                                       Int_t vtxbin) {
533   //Get background correction histogram for NSD event class
534   if(!fIsInit) {
535     AliWarning("Not initialized yet. Call Init() to remedy");
536     return 0;
537   }
538   
539   if(vtxbin > fBackground->GetNvtxBins()) {
540     AliWarning(Form("No background object for vertex bin %d", vtxbin));
541     return 0;
542   } 
543   
544   if(fBackground->GetNSDBgCorrection(det,ring,vtxbin))
545     return fBackground->GetNSDBgCorrection(det,ring,vtxbin);
546   else 
547     AliWarning("No NSD background map. You get usual one. Difference is probablyu negligible");
548   
549   return fBackground->GetBgCorrection(det,ring,vtxbin); 
550 }
551 //____________________________________________________________________
552
553 TH1F* AliFMDAnaParameters::GetDoubleHitCorrection(Int_t det, 
554                                                   Char_t ring) {
555   //Get correction for several hits in strips for p+p
556   if(!fIsInit) {
557     AliWarning("Not initialized yet. Call Init() to remedy");
558     return 0;
559   }
560   
561   return fBackground->GetDoubleHitCorrection(det,ring);
562 }
563 //_____________________________________________________________________
564 TH1F* AliFMDAnaParameters::GetSPDDeadCorrection(Int_t vtxbin) {
565                                                 
566   //Get correction for several hits in strips for p+p
567   if(!fIsInit) {
568     AliWarning("Not initialized yet. Call Init() to remedy");
569     return 0;
570   }
571   
572   return fBackground->GetSPDDeadCorrection(vtxbin);
573 }
574 //_____________________________________________________________________
575 Float_t AliFMDAnaParameters::GetEventSelectionEfficiency(Int_t vtxbin) {
576   //Get event selection efficiency object
577   if(!fIsInit) {
578     AliWarning("Not initialized yet. Call Init() to remedy");
579     return 0;
580   }
581   return fEventSelectionEfficiency->GetCorrection(vtxbin);
582
583 }
584 //_____________________________________________________________________
585 Float_t AliFMDAnaParameters::GetVtxSelectionEffFromMC() {
586
587    if(!fIsInit) {
588     AliWarning("Not initialized yet. Call Init() to remedy");
589     return 0;
590    }
591    return fEventSelectionEfficiency->GetVtxToTriggerRatio();
592
593
594 }
595 //_____________________________________________________________________
596 TH2F* AliFMDAnaParameters::GetEventSelectionEfficiency(Char_t* trig, Int_t vtxbin, Char_t ring) {
597   //Get event selection efficiency object
598   
599   TString test = trig;
600   if(!test.Contains("NSD") && !test.Contains("INEL")) {
601     AliWarning("Event selection efficiency only available for INEL and NSD");
602     return 0;
603   }
604   if(!fIsInit) {
605     AliWarning("Not initialized yet. Call Init() to remedy");
606     return 0;
607   }
608   return fEventSelectionEfficiency->GetCorrection(trig,vtxbin,ring);
609
610 }
611 //_____________________________________________________________________
612 TH1F* AliFMDAnaParameters::GetSharingEfficiency(Int_t det, Char_t ring, Int_t vtxbin) {
613   //Get sharing efficiency object
614   if(!fIsInit) {
615     AliWarning("Not initialized yet. Call Init() to remedy");
616     return 0;
617   }
618   
619   return fSharingEfficiency->GetSharingEff(det,ring,vtxbin);
620
621 }
622 //_____________________________________________________________________
623 TH1F* AliFMDAnaParameters::GetSharingEfficiencyTrVtx(Int_t det, Char_t ring, Int_t vtxbin) {
624   //Get sharing efficiency object TrVtx
625   if(!fIsInit) {
626     AliWarning("Not initialized yet. Call Init() to remedy");
627     return 0;
628   }
629   
630   return fSharingEfficiency->GetSharingEffTrVtx(det,ring,vtxbin);
631
632 }
633 //_____________________________________________________________________
634 Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) const {
635   //Get max R of ring
636   Float_t radius = 0;
637   if(ring == 'I')
638     radius = 17.2;
639   else if(ring == 'O')
640     radius = 28.0;
641   else
642     AliWarning("Unknown ring - must be I or O!");
643   
644   return radius;
645 }
646 //_____________________________________________________________________
647 Float_t AliFMDAnaParameters::GetMinR(Char_t ring) const{
648   //Get min R of ring
649   Float_t radius = 0;
650   if(ring == 'I')
651     radius = 4.5213;
652   else if(ring == 'O')
653     radius = 15.4;
654   else
655     AliWarning("Unknown ring - must be I or O!");
656   
657   return radius;
658
659 }
660 //_____________________________________________________________________
661 void AliFMDAnaParameters::SetCorners(Char_t ring) {
662   //Set corners (taken from nominal geometry)
663   if(ring == 'I') {
664     fCorner1.Set(4.9895, 15.3560);
665     fCorner2.Set(1.8007, 17.2000);
666   }
667   else {
668     fCorner1.Set(4.2231, 26.6638);
669     fCorner2.Set(1.8357, 27.9500);
670   }
671   
672 }
673 //_____________________________________________________________________
674 Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) const
675 {
676   //Get phi from sector
677   Int_t nsec = (ring == 'I' ? 20 : 40);
678   Float_t basephi = 0;
679   if(det == 1) 
680     basephi = 1.72787594; 
681   if(det == 2 && ring == 'I')
682     basephi = 0.15707963;
683   if(det == 2 && ring == 'O')
684     basephi = 0.078539818;
685   if(det == 3 && ring == 'I')
686     basephi = 2.984513044;
687   if(det == 3 && ring == 'O')
688     basephi = 3.06305289;
689   
690   Float_t step = 2*TMath::Pi() / nsec;
691   Float_t phi = 0;
692   if(det == 3)
693     phi = basephi - sec*step;
694   else
695     phi = basephi + sec*step;
696   
697   if(phi < 0) 
698     phi = phi +2*TMath::Pi();
699   if(phi > 2*TMath::Pi() )
700     phi = phi - 2*TMath::Pi();
701   
702   return phi;
703 }
704 //_____________________________________________________________________
705 Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) const
706 {
707   //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta)
708   Float_t   rad       = GetMaxR(ring)-GetMinR(ring);
709   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
710   Float_t   segment   = rad / nStrips;
711   Float_t   r         = GetMinR(ring) + segment*strip;
712   Float_t   z         = 0;
713   Int_t hybrid = sec / 2;
714   
715   if(det == 1) {
716     if(!(hybrid%2)) z = 320.266; else z = 319.766;
717   }
718   if(det == 2 && ring == 'I' ) {
719     if(!(hybrid%2)) z = 83.666; else z = 83.166;
720   }
721   if(det == 2 && ring == 'O' ) {
722     if(!(hybrid%2)) z = 74.966; else z = 75.466;
723   }
724   if(det == 3 && ring == 'I' ) {
725     if(!(hybrid%2)) z = -63.066; else z = -62.566;
726   }
727   if(det == 3 && ring == 'O' ) {
728     if(!(hybrid%2)) z = -74.966; else z = -75.466;
729   }
730   
731   //std::cout<<det<<"   "<<ring<<"   "<<sec<<"   "<<hybrid<<"    "<<z<<std::endl;
732   
733   // Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
734   Float_t   theta = TMath::ATan2(r,z-zvtx);
735   Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
736   
737   return eta;
738 }
739
740 //_____________________________________________________________________
741
742 Bool_t AliFMDAnaParameters::GetVertex(const AliESDEvent* esd, Double_t* vertexXYZ) 
743 {
744   //Get the vertex from the ESD
745   const AliESDVertex* vertex = 0;
746   vertex = esd->GetPrimaryVertexSPD();
747   
748   if(vertex)
749     vertex->GetXYZ(vertexXYZ);
750
751   //if(vertexXYZ[0] == 0 || vertexXYZ[1] == 0 )
752   //  return kFALSE;
753  
754   if(vertex->GetNContributors() <= 0)
755     return kFALSE;
756   
757   if(vertex->GetZRes() > 0.1 ) 
758     return kFALSE;
759   
760   return vertex->GetStatus();
761   
762 }
763 //____________________________________________________________________
764 void AliFMDAnaParameters::SetTriggerStatus(const AliESDEvent *esd) {
765
766   //ULong64_t triggerMask = esd->GetTriggerMask();
767   
768   AliPhysicsSelection* centralPhysicsSelection = (AliPhysicsSelection*)((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetEventSelection();
769   
770   if(!centralPhysicsSelection && !fPhysicsSelection) {
771     
772     std::cout<<"Creating AliPhysicsSelection object due to absence of central object"<<std::endl;
773     fPhysicsSelection = new AliPhysicsSelection;
774     fPhysicsSelection->SetAnalyzeMC(!fRealData);
775     // fPhysicsSelection->SetUseBXNumbers(kFALSE);
776     
777     
778     AliBackgroundSelection* backgroundSelection = new AliBackgroundSelection("bg","bg");
779     backgroundSelection->Init();
780     fPhysicsSelection->AddBackgroundIdentification(backgroundSelection);
781     
782   }
783   
784   TString triggers = esd->GetFiredTriggerClasses();
785   
786   AliTriggerAnalysis tAna;
787   
788   fTriggerInel  = kFALSE;
789   fTriggerNSD   = kFALSE;
790   fTriggerEmpty = kFALSE;
791   
792   UInt_t inel = kFALSE;
793   
794   if(centralPhysicsSelection) {
795     inel= ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
796   }
797   else
798     inel = fPhysicsSelection->IsCollisionCandidate(esd);
799   
800   
801   
802   if(fInelGtZero && inel) {
803     const AliMultiplicity* spdmult = esd->GetMultiplicity();
804     Int_t nCentralTracklets = 0;
805     Int_t j = 0;
806     while( nCentralTracklets < 1 && j< spdmult->GetNumberOfTracklets() ) {
807       if(TMath::Abs(spdmult->GetEta(j)) < 1) nCentralTracklets++;
808       j++;
809     }
810     if(nCentralTracklets < 1) inel = kFALSE;      
811   }
812   
813   //std::cout<<fTriggerInel<<std::endl;
814   if(inel) {
815     fTriggerInel = kTRUE;
816   }
817   
818   Bool_t nsd = kFALSE;
819   if(fUseBuiltInNSD) { 
820     if ((tAna.FMDTrigger(esd, AliTriggerAnalysis::kASide) || tAna.V0Trigger(esd, AliTriggerAnalysis::kASide, kFALSE) == AliTriggerAnalysis::kV0BB) && (tAna.FMDTrigger(esd, AliTriggerAnalysis::kCSide) || tAna.V0Trigger(esd, AliTriggerAnalysis::kCSide, kFALSE) == AliTriggerAnalysis::kV0BB))
821       nsd = kTRUE;
822   }
823   else nsd = tAna.IsOfflineTriggerFired(esd,AliTriggerAnalysis::kNSD1);
824   
825   if(fTriggerInel && nsd) {
826     fTriggerNSD = kTRUE;
827   }
828   if(triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
829     fTriggerEmpty = kTRUE;
830   }
831   
832   
833   /*switch (fTrigger) {
834   case kMB1: {
835     if( fPhysicsSelection->IsCollisionCandidate(esd)) {
836       fTriggerInel = kTRUE;
837     }
838     break;
839     
840   }
841   case kMB2: { 
842     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
843       return kTRUE;
844     break;
845   }
846   case kSPDFASTOR: {
847     if (triggerMask & spdFO)
848       return kTRUE;
849     break;
850   }
851   case kNOCTP: {
852     return kTRUE;
853     break;
854   }
855   case kEMPTY: {
856   if(triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) 
857   return kTRUE;
858     break;
859   }
860   case kNSD: {
861     if(fPhysicsSelection->IsCollisionCandidate(esd) && tAna.IsOfflineTriggerFired(esd,AliTriggerAnalysis::kNSD1))
862       return kTRUE;
863     break;
864   }
865     
866   }//switch
867   */
868   
869 }
870 /*
871 //____________________________________________________________________
872 Bool_t AliFMDAnaParameters::IsEventTriggered(const AliESDEvent *esd, Trigger trig) {
873   //Did we have trig trigger ?
874   Trigger old = fTrigger;
875   fTrigger = trig;
876   Bool_t retval = IsEventTriggered(esd);
877   fTrigger = old;
878   return retval;
879
880 }
881 */
882 //____________________________________________________________________
883 Bool_t AliFMDAnaParameters::IsEventTriggered(Trigger trigger) {
884   // check if the event was triggered
885   
886   if (fCentralSelection) return kTRUE;
887   switch (trigger) {
888   
889   case kMB1:
890     return fTriggerInel;
891     break;
892   case kNSD:
893     return fTriggerNSD;
894     break;
895   case kEMPTY:
896     return fTriggerEmpty;
897     break;
898   case kNOCTP:
899     return kTRUE;
900     break;
901   default:
902     AliWarning("Trigger not implemented!!!");
903     break;
904     
905     
906   }
907   return kFALSE;
908   
909 }
910     
911 //____________________________________________________________________
912 Float_t 
913 AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)  
914 {
915   //Get length of a strip
916   
917   Float_t rad        = GetMaxR(ring)-GetMinR(ring);
918   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
919   Float_t segment    = rad / nStrips;
920   
921   //TVector2* corner1  = fmdring.GetVertex(2);  
922   // TVector2* corner2  = fmdring.GetVertex(3);
923   
924   SetCorners(ring);
925   /*
926   std::cout<<GetMaxR(ring)<<"   "<<fmdring.GetMaxR()<<std::endl;
927   std::cout<<GetMinR(ring)<<"   "<<fmdring.GetMinR()<<std::endl;
928   std::cout<<corner1->X()<<"   "<<fCorner1.X()<<std::endl;
929   std::cout<<corner2->X()<<"   "<<fCorner2.X()<<std::endl;
930   std::cout<<corner1->Y()<<"   "<<fCorner1.Y()<<std::endl;
931   std::cout<<corner2->Y()<<"   "<<fCorner2.Y()<<std::endl;*/
932   Float_t slope      = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
933   Float_t constant   = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
934   Float_t radius     = GetMinR(ring) + strip*segment;
935   
936   Float_t d          = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
937   
938   Float_t arclength  = GetBaseStripLength(ring,strip);
939   if(d>0) {
940     
941     Float_t x        = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
942     Float_t y        = slope*x + constant;
943     Float_t theta    = TMath::ATan2(x,y);
944     
945     if(x < fCorner1.X() && y > fCorner1.Y()) {
946       arclength = radius*theta;                        //One sector since theta is by definition half-hybrid
947       
948     }
949     
950   }
951   
952   return arclength;
953   
954   
955 }
956 //____________________________________________________________________
957 Float_t 
958 AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)  
959 {  
960   //Get length of strip assuming that corners are not cut away
961   Float_t rad             = GetMaxR(ring)-GetMinR(ring);
962   Float_t nStrips         = (ring == 'I' ? 512 : 256);
963   Float_t nSec            = (ring == 'I' ? 20 : 40);
964   Float_t segment         = rad / nStrips;
965   Float_t basearc         = 2*TMath::Pi() / (0.5*nSec); // One hybrid: 36 degrees inner, 18 outer
966   Float_t radius          = GetMinR(ring) + strip*segment;
967   Float_t basearclength   = 0.5*basearc * radius;                // One sector   
968   
969   return basearclength;
970 }
971 //____________________________________________________________________
972 Int_t    AliFMDAnaParameters::GetFirstEtaBinFromMap(Int_t vtxbin, Int_t det, Char_t ring)
973 {
974   TH2F* hBg = GetBackgroundCorrection(det,ring,vtxbin);
975   
976   if(det == 0) return hBg->GetXaxis()->FindBin(-1.95);
977   
978   Int_t firstbin = -1;
979   Int_t nNonZeroFirst = 0;
980   
981   for(Int_t i=1;i<=hBg->GetNbinsX();i++) {
982     if(nNonZeroFirst == fNumberOfEtaBinsToCut && firstbin==-1) firstbin = i;
983     
984     for(Int_t j=1;j<=hBg->GetNbinsY();j++) {
985       
986       Float_t value = hBg->GetBinContent(i,j);
987       
988       if(value > 0.001 && nNonZeroFirst<fNumberOfEtaBinsToCut)
989         {nNonZeroFirst++; break;}
990       
991       
992     }
993   }
994   
995   return firstbin;
996
997 }
998 //____________________________________________________________________
999 Int_t    AliFMDAnaParameters::GetLastEtaBinFromMap(Int_t vtxbin, Int_t det, Char_t ring)
1000 {
1001   TH2F* hBg = GetBackgroundCorrection(det,ring,vtxbin);
1002   Int_t lastbin=-1;
1003   Int_t nNonZeroLast = 0;
1004   
1005   if(det == 0) return hBg->GetXaxis()->FindBin(1.95);
1006   
1007   for(Int_t i=hBg->GetNbinsX();i>0;i--) {
1008     if(nNonZeroLast == fNumberOfEtaBinsToCut && lastbin==-1) lastbin = i;
1009     
1010     for(Int_t j=1;j<=hBg->GetNbinsY();j++) {
1011       
1012       Float_t value = hBg->GetBinContent(i,j);
1013       
1014       if(value > 0.001 && nNonZeroLast<fNumberOfEtaBinsToCut)
1015         {nNonZeroLast++; break; }
1016       
1017       
1018     }
1019   }
1020   
1021   return lastbin;
1022 }
1023
1024 //____________________________________________________________________
1025 Int_t    AliFMDAnaParameters::GetFirstEtaBinToInclude(Int_t vtxbin, Int_t det, Char_t ring)
1026 {
1027   Int_t ringNumber = (ring == 'I' ? 0 : 1);
1028   return fEtaLowBinLimits.GetBinContent(det,ringNumber,vtxbin);
1029
1030 }
1031
1032 //____________________________________________________________________
1033 Int_t    AliFMDAnaParameters::GetLastEtaBinToInclude(Int_t vtxbin, Int_t det, Char_t ring)
1034 {
1035   Int_t ringNumber = (ring == 'I' ? 0 : 1);
1036   return fEtaHighBinLimits.GetBinContent(det,ringNumber,vtxbin);
1037   
1038 }
1039 //____________________________________________________________________
1040 //
1041 // EOF
1042 //