PWG2/FORWARD/analysis/AliFMDAnaParameters.cxx

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