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