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