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