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