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   TString s(sys);
 304   s.ToLower();
 305   if      (s.Contains("p-p")   || s.Contains("pp"))   fSpecies = kPP;
 306   else if (s.Contains("pb-pb") || s.Contains("pbpb")) fSpecies = kPbPb;
 307   else if (s.Contains("a-a")   || s.Contains("aa"))   fSpecies = kPbPb;
 308 }
 309
 310 //____________________________________________________________________
 311 void AliFMDAnaParameters::SetParametersFromESD(AliESDEvent* esd)
 312 {
 313
 314
 315   SetCollisionSystem(esd->GetBeamType());
 316
 317   Float_t energy   = esd->GetBeamEnergy();
 318   // Correct to center of mass per nucleon (cmsNN) - LHC gives it as
 319   // cmsNN * Z / 2
 320   if (fSpecies == kPbPb) energy = energy / 208 * 82;
 321   SetEnergy(2*energy);
 322   SetMagField(esd->GetMagneticField());
 323
 324   // Float_t magfield = esd->GetCurrentL3();
 325   // if (TMath::Abs(magfield - 30000.) < 10 ) fMagField = k5G;
 326   // if (TMath::Abs(magfield + 30000.) < 10 ) fMagField = k5Gnegative;
 327   // if (TMath::Abs(magfield) < 10 )          fMagField = k0G;
 328
 329
 330
 331   Init(kTRUE);
 332
 333 }
 334
 335 //____________________________________________________________________
 336
 337 void AliFMDAnaParameters::PrintStatus(Bool_t showpaths) const
 338 {
 339   //Print current status
 340   TString energystring;
 341   switch(fEnergy) {
 342   case k900:
 343     energystring.Form("900 GeV");   break;
 344   case k2400:
 345     energystring.Form("2400 GeV"); break;
 346   case k2750:
 347     energystring.Form("2750 GeV"); break;
 348   case k5500:
 349     energystring.Form("5500 GeV"); break;
 350   case k7000:
 351     energystring.Form("7000 GeV");  break;
 352   case k10000:
 353     energystring.Form("10000 GeV"); break;
 354   case k14000:
 355     energystring.Form("14000 GeV"); break;
 356   default:
 357     energystring.Form("invalid energy"); break;
 358   }
 359   TString triggerstring;
 360   switch(fTrigger) {
 361   case kMB1:
 362     triggerstring.Form("Minimum bias 1");   break;
 363   case kMB2:
 364     triggerstring.Form("Minimum bias 2");   break;
 365   case kSPDFASTOR:
 366     triggerstring.Form("SPD FAST OR");   break;
 367   case kNOCTP:
 368     triggerstring.Form("NO TRIGGER TEST");   break;
 369   default:
 370     energystring.Form("invalid trigger"); break;
 371   }
 372   TString magstring;
 373   switch(fMagField) {
 374   case k5G:
 375     magstring.Form("+5 kGaus");   break;
 376   case k0G:
 377     magstring.Form("0 kGaus");   break;
 378   case k5Gnegative:
 379     magstring.Form("-5 kGaus");   break;
 380   default:
 381     magstring.Form("invalid mag field %d", fMagField); break;
 382   }
 383   TString collsystemstring;
 384   switch(fSpecies) {
 385   case kPP:
 386     collsystemstring.Form("p-p");   break;
 387   case kPbPb:
 388     collsystemstring.Form("Pb-Pb");   break;
 389   default:
 390     collsystemstring.Form("invalid collision system");   break;
 391   }
 392
 393   TString datastring;
 394   switch(fRealData) {
 395   case kTRUE:
 396     datastring.Form("Nature"); break;
 397   case kFALSE:
 398     datastring.Form("MC"); break;
 399   default:
 400     datastring.Form("Unknown"); break ;
 401   }
 402
 403   TString InelString;
 404   if(fInelGtZero) InelString = "INEL > 0";
 405   else InelString = "INEL";
 406
 407   std::cout<<"Energy       = "<<energystring.Data()<<std::endl;
 408   std::cout<<"Trigger      = "<<triggerstring.Data()<<std::endl;
 409   std::cout<<"Mag Field    = "<<magstring.Data()<<std::endl;
 410   std::cout<<"Coll System  = "<<collsystemstring.Data()<<std::endl;
 411   std::cout<<"Data origin  = "<<datastring.Data()<<std::endl;
 412   std::cout<<"Basic trigger: "<<InelString.Data()<<std::endl;
 413
 414   if (showpaths) {
 415     TString bg = GetPath(fgkBackgroundID);
 416     TString es = GetPath(fgkEventSelectionEffID);
 417     TString ed = GetPath(fgkEnergyDistributionID);
 418     TString me = GetPath(fgkSharingEffID);
 419     std::cout << "2nd maps:     " << bg << "\n"
 420               << "Event sel.:   " << es << "\n"
 421               << "Energy dist.: " << ed << "\n"
 422               << "Merge eff.:   " << me << std::endl;
 423   }
 424
 425 }
 426
 427 //____________________________________________________________________
 428 Float_t AliFMDAnaParameters::GetVtxCutZ() {
 429   //Get the z vtx cut in analysis
 430   if(!fIsInit) {
 431     AliWarning("Not initialized yet. Call Init() to remedy");
 432     return -1;
 433   }
 434
 435   return fBackground->GetVtxCutZ();
 436 }
 437
 438 //____________________________________________________________________
 439 Int_t AliFMDAnaParameters::GetNvtxBins() {
 440   //Get number of vtx bins
 441   if(!fIsInit) {
 442     AliWarning("Not initialized yet. Call Init() to remedy");
 443     return -1;
 444   }
 445
 446   return fBackground->GetNvtxBins();
 447 }
 448 //____________________________________________________________________
 449 TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
 450
 451   return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
 452 }
 453 //____________________________________________________________________
 454 TH1F* AliFMDAnaParameters::GetEmptyEnergyDistribution(Int_t det, Char_t ring) {
 455
 456   return fEnergyDistribution->GetEmptyEnergyDistribution(det, ring);
 457 }
 458 //____________________________________________________________________
 459 TH1F* AliFMDAnaParameters::GetRingEnergyDistribution(Int_t det, Char_t ring) {
 460
 461   return fEnergyDistribution->GetRingEnergyDistribution(det, ring);
 462 }
 463 //____________________________________________________________________
 464 Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
 465   //Get sigma of Landau fits
 466   if(!fIsInit) {
 467     AliWarning("Not initialized yet. Call Init() to remedy");
 468     return 0;
 469   }
 470
 471   TH1F* hEnergyDist       = GetEnergyDistribution(det,ring, eta);
 472   TF1*  fitFunc           = hEnergyDist->GetFunction("FMDfitFunc");
 473   if(!fitFunc) {
 474     //AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
 475     return 1024;
 476   }
 477   Float_t sigma           = fitFunc->GetParameter(2);
 478   return sigma;
 479 }
 480
 481
 482 //____________________________________________________________________
 483 Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
 484   //Get MPV of landau fits
 485   if(!fIsInit) {
 486     AliWarning("Not initialized yet. Call Init() to remedy");
 487     return 0;
 488   }
 489   //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
 490   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 491   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 492
 493   if(!fitFunc) {
 494     AliWarning(Form("No function for FMD%d%c, eta %f (%d)",det,ring,eta,
 495                     GetEtaBin(eta)));
 496     return 1024;
 497   }
 498
 499   Float_t mpv           = fitFunc->GetParameter(1);
 500   return mpv;
 501 }
 502 //____________________________________________________________________
 503 Float_t AliFMDAnaParameters::GetConstant(Int_t det, Char_t ring, Float_t eta) {
 504   //Get constant parameter of Landau fits
 505   if(!fIsInit) {
 506     AliWarning("Not initialized yet. Call Init() to remedy");
 507     return 0;
 508   }
 509
 510   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 511   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 512   if(!fitFunc) {
 513     AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
 514     return 0;
 515   }
 516
 517   Float_t mpv           = fitFunc->GetParameter(0);
 518   return mpv;
 519 }
 520 //____________________________________________________________________
 521 Float_t
 522 AliFMDAnaParameters::Get2MIPWeight(Int_t det, Char_t ring, Float_t eta)
 523 {
 524   //Get 2 MIP weights of convoluted Landau fits
 525   if(!fIsInit) {
 526     AliWarning("Not initialized yet. Call Init() to remedy");
 527     return 0;
 528   }
 529
 530   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 531   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 532   if(!fitFunc) return 0;
 533   Float_t twoMIPweight    = fitFunc->GetParameter(3);
 534
 535
 536
 537   if(twoMIPweight < 1e-05)
 538     twoMIPweight = 0;
 539
 540   return twoMIPweight;
 541 }
 542 //____________________________________________________________________
 543 Float_t
 544 AliFMDAnaParameters::Get3MIPWeight(Int_t det, Char_t ring, Float_t eta)
 545 {
 546   //Get 3 MIP weights of convoluted Landau fits
 547   if(!fIsInit) {
 548     AliWarning("Not initialized yet. Call Init() to remedy");
 549     return 0;
 550   }
 551
 552   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 553   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 554   if(!fitFunc) return 0;
 555   Float_t threeMIPweight    = fitFunc->GetParameter(4);
 556
 557   if(threeMIPweight < 1e-05)
 558     threeMIPweight = 0;
 559
 560   Float_t twoMIPweight    = fitFunc->GetParameter(3);
 561
 562   if(twoMIPweight < 1e-05)
 563     threeMIPweight = 0;
 564
 565   return threeMIPweight;
 566 }
 567 //____________________________________________________________________
 568 Int_t AliFMDAnaParameters::GetNetaBins()
 569 {
 570   return GetBackgroundCorrection(1,'I',5)->GetNbinsX();
 571 }
 572 //____________________________________________________________________
 573 Float_t AliFMDAnaParameters::GetEtaMin()
 574 {
 575   return GetBackgroundCorrection(1,'I',5)->GetXaxis()->GetXmin();
 576 }
 577 //____________________________________________________________________
 578 Float_t AliFMDAnaParameters::GetEtaMax()
 579 {
 580   return GetBackgroundCorrection(1,'I',5)->GetXaxis()->GetXmax();
 581 }
 582 //____________________________________________________________________
 583 Int_t AliFMDAnaParameters::GetEtaBin(Float_t eta)
 584 {
 585   TAxis testaxis(GetNetaBins(),GetEtaMin(),GetEtaMax());
 586   Int_t binnumber = testaxis.FindBin(eta) ;
 587
 588   return binnumber;
 589 }
 590 //____________________________________________________________________
 591
 592 TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det,
 593                                                    Char_t ring,
 594                                                    Int_t vtxbin) {
 595   //Get background correction histogram
 596   if(!fIsInit) {
 597     AliWarning("Not initialized yet. Call Init() to remedy");
 598     return 0;
 599   }
 600
 601
 602
 603   if(vtxbin > fBackground->GetNvtxBins()) {
 604     AliWarning(Form("No background object for vertex bin %d", vtxbin));
 605     return 0;
 606   }
 607
 608   return fBackground->GetBgCorrection(det,ring,vtxbin);
 609 }
 610 //____________________________________________________________________
 611 TH2F* AliFMDAnaParameters::GetBackgroundCorrectionNSD(Int_t det,
 612                                                       Char_t ring,
 613                                                       Int_t vtxbin) {
 614   //Get background correction histogram for NSD event class
 615   if(!fIsInit) {
 616     AliWarning("Not initialized yet. Call Init() to remedy");
 617     return 0;
 618   }
 619
 620   if(vtxbin > fBackground->GetNvtxBins()) {
 621     AliWarning(Form("No background object for vertex bin %d", vtxbin));
 622     return 0;
 623   }
 624
 625   if(fBackground->GetNSDBgCorrection(det,ring,vtxbin))
 626     return fBackground->GetNSDBgCorrection(det,ring,vtxbin);
 627   else
 628     AliWarning("No NSD background map. You get usual one. "
 629                "Difference is probably negligible");
 630
 631   return fBackground->GetBgCorrection(det,ring,vtxbin);
 632 }
 633 //____________________________________________________________________
 634
 635 TH1F* AliFMDAnaParameters::GetDoubleHitCorrection(Int_t det,
 636                                                   Char_t ring) {
 637   //Get correction for several hits in strips for p+p
 638   if(!fIsInit) {
 639     AliWarning("Not initialized yet. Call Init() to remedy");
 640     return 0;
 641   }
 642
 643   return fBackground->GetDoubleHitCorrection(det,ring);
 644 }
 645 //_____________________________________________________________________
 646 TH1F* AliFMDAnaParameters::GetSPDDeadCorrection(Int_t vtxbin) {
 647
 648   //Get correction for several hits in strips for p+p
 649   if(!fIsInit) {
 650     AliWarning("Not initialized yet. Call Init() to remedy");
 651     return 0;
 652   }
 653
 654   return fBackground->GetSPDDeadCorrection(vtxbin);
 655 }
 656 //_____________________________________________________________________
 657 TH1F* AliFMDAnaParameters::GetFMDDeadCorrection(Int_t vtxbin) {
 658
 659   //Get correction for several hits in strips for p+p
 660   if(!fIsInit) {
 661     AliWarning("Not initialized yet. Call Init() to remedy");
 662     return 0;
 663   }
 664
 665   return fBackground->GetFMDDeadCorrection(vtxbin);
 666 }
 667 //_____________________________________________________________________
 668 Float_t AliFMDAnaParameters::GetEventSelectionEfficiency(Int_t vtxbin) {
 669   //Get event selection efficiency object
 670   if(!fIsInit) {
 671     AliWarning("Not initialized yet. Call Init() to remedy");
 672     return 0;
 673   }
 674   return fEventSelectionEfficiency->GetCorrection(vtxbin);
 675
 676 }
 677 //_____________________________________________________________________
 678 Float_t AliFMDAnaParameters::GetVtxSelectionEffFromMC() {
 679
 680    if(!fIsInit) {
 681     AliWarning("Not initialized yet. Call Init() to remedy");
 682     return 0;
 683    }
 684    return fEventSelectionEfficiency->GetVtxToTriggerRatio();
 685
 686
 687 }
 688 //_____________________________________________________________________
 689 TH2F* AliFMDAnaParameters::GetEventSelectionEfficiency(TString trig, Int_t vtxbin, Char_t ring) {
 690   //Get event selection efficiency object
 691
 692   //TString test = trig;
 693   if(!trig.Contains("NSD") && !trig.Contains("INEL")) {
 694     AliWarning("Event selection efficiency only available for INEL and NSD");
 695     return 0;
 696   }
 697   if(!fIsInit) {
 698     AliWarning("Not initialized yet. Call Init() to remedy");
 699     return 0;
 700   }
 701   return fEventSelectionEfficiency->GetCorrection(trig,vtxbin,ring);
 702
 703 }
 704 //_____________________________________________________________________
 705 TH1F* AliFMDAnaParameters::GetSharingEfficiency(Int_t det, Char_t ring, Int_t vtxbin) {
 706   //Get sharing efficiency object
 707   if(!fIsInit) {
 708     AliWarning("Not initialized yet. Call Init() to remedy");
 709     return 0;
 710   }
 711
 712   return fSharingEfficiency->GetSharingEff(det,ring,vtxbin);
 713
 714 }
 715 //_____________________________________________________________________
 716 TH1F* AliFMDAnaParameters::GetSharingEfficiencyTrVtx(Int_t det, Char_t ring, Int_t vtxbin) {
 717   //Get sharing efficiency object TrVtx
 718   if(!fIsInit) {
 719     AliWarning("Not initialized yet. Call Init() to remedy");
 720     return 0;
 721   }
 722
 723   return fSharingEfficiency->GetSharingEffTrVtx(det,ring,vtxbin);
 724
 725 }
 726 //_____________________________________________________________________
 727 Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) const {
 728   //Get max R of ring
 729   Float_t radius = 0;
 730   if(ring == 'I')
 731     radius = 17.2;
 732   else if(ring == 'O')
 733     radius = 28.0;
 734   else
 735     AliWarning("Unknown ring - must be I or O!");
 736
 737   return radius;
 738 }
 739 //_____________________________________________________________________
 740 Float_t AliFMDAnaParameters::GetMinR(Char_t ring) const{
 741   //Get min R of ring
 742   Float_t radius = 0;
 743   if(ring == 'I')
 744     radius = 4.5213;
 745   else if(ring == 'O')
 746     radius = 15.4;
 747   else
 748     AliWarning("Unknown ring - must be I or O!");
 749
 750   return radius;
 751
 752 }
 753 //_____________________________________________________________________
 754 void AliFMDAnaParameters::SetCorners(Char_t ring) {
 755   //Set corners (taken from nominal geometry)
 756   if(ring == 'I') {
 757     fCorner1.Set(4.9895, 15.3560);
 758     fCorner2.Set(1.8007, 17.2000);
 759   }
 760   else {
 761     fCorner1.Set(4.2231, 26.6638);
 762     fCorner2.Set(1.8357, 27.9500);
 763   }
 764
 765 }
 766 //_____________________________________________________________________
 767 Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) const
 768 {
 769   //Get phi from sector
 770   Int_t nsec = (ring == 'I' ? 20 : 40);
 771   Float_t basephi = 0;
 772   if(det == 1)
 773     basephi = 1.72787594;
 774   if(det == 2 && ring == 'I')
 775     basephi = 0.15707963;
 776   if(det == 2 && ring == 'O')
 777     basephi = 0.078539818;
 778   if(det == 3 && ring == 'I')
 779     basephi = 2.984513044;
 780   if(det == 3 && ring == 'O')
 781     basephi = 3.06305289;
 782
 783   Float_t step = 2*TMath::Pi() / nsec;
 784   Float_t phi = 0;
 785   if(det == 3)
 786     phi = basephi - sec*step;
 787   else
 788     phi = basephi + sec*step;
 789
 790   if(phi < 0)
 791     phi = phi +2*TMath::Pi();
 792   if(phi > 2*TMath::Pi() )
 793     phi = phi - 2*TMath::Pi();
 794
 795   return phi;
 796 }
 797 //_____________________________________________________________________
 798 Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) const
 799 {
 800   //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta)
 801   Float_t   rad       = GetMaxR(ring)-GetMinR(ring);
 802   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
 803   Float_t   segment   = rad / nStrips;
 804   Float_t   r         = GetMinR(ring) + segment*strip;
 805   Float_t   z         = 0;
 806   Int_t hybrid = sec / 2;
 807
 808   if(det == 1) {
 809     if(!(hybrid%2)) z = 320.266; else z = 319.766;
 810   }
 811   if(det == 2 && ring == 'I' ) {
 812     if(!(hybrid%2)) z = 83.666; else z = 83.166;
 813   }
 814   if(det == 2 && ring == 'O' ) {
 815     if(!(hybrid%2)) z = 74.966; else z = 75.466;
 816   }
 817   if(det == 3 && ring == 'I' ) {
 818     if(!(hybrid%2)) z = -63.066; else z = -62.566;
 819   }
 820   if(det == 3 && ring == 'O' ) {
 821     if(!(hybrid%2)) z = -74.966; else z = -75.466;
 822   }
 823
 824   //std::cout<<det<<"   "<<ring<<"   "<<sec<<"   "<<hybrid<<"    "<<z<<std::endl;
 825
 826   // Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 827   Float_t   theta = TMath::ATan2(r,z-zvtx);
 828   Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 829
 830   return eta;
 831 }
 832
 833 //_____________________________________________________________________
 834
 835 Bool_t AliFMDAnaParameters::GetVertex(const AliESDEvent* esd, Double_t* vertexXYZ)
 836 {
 837   //Get the vertex from the ESD
 838   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
 839
 840   if (!vertex) return kFALSE;
 841
 842   vertex->GetXYZ(vertexXYZ);
 843
 844   //if(vertexXYZ[0] == 0 || vertexXYZ[1] == 0 )
 845   //  return kFALSE;
 846
 847   if(vertex->GetNContributors() <= 0)
 848     return kFALSE;
 849
 850   if(vertex->GetZRes() > 0.1 )
 851     return kFALSE;
 852
 853   return vertex->GetStatus();
 854
 855 }
 856 //____________________________________________________________________
 857 void AliFMDAnaParameters::SetTriggerStatus(const AliESDEvent *esd) {
 858
 859   //ULong64_t triggerMask = esd->GetTriggerMask();
 860
 861   AliPhysicsSelection* centralPhysicsSelection = (AliPhysicsSelection*)((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetEventSelection();
 862
 863   if(!centralPhysicsSelection && !fPhysicsSelection) {
 864
 865     std::cout<<"Creating AliPhysicsSelection object due to absence of central object"<<std::endl;
 866     fPhysicsSelection = new AliPhysicsSelection;
 867     fPhysicsSelection->SetAnalyzeMC(!fRealData);
 868     // fPhysicsSelection->SetUseBXNumbers(kFALSE);
 869
 870
 871     AliBackgroundSelection* backgroundSelection = new AliBackgroundSelection("bg","bg");
 872     backgroundSelection->Init();
 873     fPhysicsSelection->AddBackgroundIdentification(backgroundSelection);
 874
 875   }
 876
 877   TString triggers = esd->GetFiredTriggerClasses();
 878
 879   AliTriggerAnalysis tAna;
 880
 881   fTriggerInel  = kFALSE;
 882   fTriggerNSD   = kFALSE;
 883   fTriggerEmpty = kFALSE;
 884
 885   UInt_t inel = kFALSE;
 886
 887   if(centralPhysicsSelection) {
 888     inel= ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 889   }
 890   else
 891     inel = fPhysicsSelection->IsCollisionCandidate(esd);
 892
 893
 894
 895   if(fInelGtZero && inel) {
 896     const AliMultiplicity* spdmult = esd->GetMultiplicity();
 897     Int_t nCentralTracklets = 0;
 898     Int_t j = 0;
 899     while( nCentralTracklets < 1 && j< spdmult->GetNumberOfTracklets() ) {
 900       if(TMath::Abs(spdmult->GetEta(j)) < 1) nCentralTracklets++;
 901       j++;
 902     }
 903     if(nCentralTracklets < 1) inel = kFALSE;
 904   }
 905
 906   //std::cout<<fTriggerInel<<std::endl;
 907   if(inel) {
 908     fTriggerInel = kTRUE;
 909   }
 910
 911   Bool_t nsd = kFALSE;
 912   if(fUseBuiltInNSD) {
 913     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))
 914       nsd = kTRUE;
 915   }
 916   else nsd = tAna.IsOfflineTriggerFired(esd,AliTriggerAnalysis::kNSD1);
 917
 918   if(fTriggerInel && nsd) {
 919     fTriggerNSD = kTRUE;
 920   }
 921   if(triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
 922     fTriggerEmpty = kTRUE;
 923   }
 924
 925
 926   /*switch (fTrigger) {
 927   case kMB1: {
 928     if( fPhysicsSelection->IsCollisionCandidate(esd)) {
 929       fTriggerInel = kTRUE;
 930     }
 931     break;
 932
 933   }
 934   case kMB2: {
 935     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 936       return kTRUE;
 937     break;
 938   }
 939   case kSPDFASTOR: {
 940     if (triggerMask & spdFO)
 941       return kTRUE;
 942     break;
 943   }
 944   case kNOCTP: {
 945     return kTRUE;
 946     break;
 947   }
 948   case kEMPTY: {
 949   if(triggers.Contains("CBEAMB-ABCE-NOPF-ALL"))
 950   return kTRUE;
 951     break;
 952   }
 953   case kNSD: {
 954     if(fPhysicsSelection->IsCollisionCandidate(esd) && tAna.IsOfflineTriggerFired(esd,AliTriggerAnalysis::kNSD1))
 955       return kTRUE;
 956     break;
 957   }
 958
 959   }//switch
 960   */
 961
 962 }
 963 /*
 964 //____________________________________________________________________
 965 Bool_t AliFMDAnaParameters::IsEventTriggered(const AliESDEvent *esd, Trigger trig) {
 966   //Did we have trig trigger ?
 967   Trigger old = fTrigger;
 968   fTrigger = trig;
 969   Bool_t retval = IsEventTriggered(esd);
 970   fTrigger = old;
 971   return retval;
 972
 973 }
 974 */
 975 //____________________________________________________________________
 976 Bool_t AliFMDAnaParameters::IsEventTriggered(Trigger trigger) {
 977   // check if the event was triggered
 978
 979   if (fCentralSelection) return kTRUE;
 980   switch (trigger) {
 981
 982   case kMB1:
 983     return fTriggerInel;
 984     break;
 985   case kNSD:
 986     return fTriggerNSD;
 987     break;
 988   case kEMPTY:
 989     return fTriggerEmpty;
 990     break;
 991   case kNOCTP:
 992     return kTRUE;
 993     break;
 994   default:
 995     AliWarning("Trigger not implemented!!!");
 996     break;
 997
 998
 999   }
1000   return kFALSE;
1001
1002 }
1003
1004 //____________________________________________________________________
1005 Float_t
1006 AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)
1007 {
1008   //Get length of a strip
1009
1010   Float_t rad        = GetMaxR(ring)-GetMinR(ring);
1011   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
1012   Float_t segment    = rad / nStrips;
1013
1014   //TVector2* corner1  = fmdring.GetVertex(2);
1015   // TVector2* corner2  = fmdring.GetVertex(3);
1016
1017   SetCorners(ring);
1018   /*
1019   std::cout<<GetMaxR(ring)<<"   "<<fmdring.GetMaxR()<<std::endl;
1020   std::cout<<GetMinR(ring)<<"   "<<fmdring.GetMinR()<<std::endl;
1021   std::cout<<corner1->X()<<"   "<<fCorner1.X()<<std::endl;
1022   std::cout<<corner2->X()<<"   "<<fCorner2.X()<<std::endl;
1023   std::cout<<corner1->Y()<<"   "<<fCorner1.Y()<<std::endl;
1024   std::cout<<corner2->Y()<<"   "<<fCorner2.Y()<<std::endl;*/
1025   Float_t slope      = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
1026   Float_t constant   = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
1027   Float_t radius     = GetMinR(ring) + strip*segment;
1028
1029   Float_t d          = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
1030
1031   Float_t arclength  = GetBaseStripLength(ring,strip);
1032   if(d>0) {
1033
1034     Float_t x        = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
1035     Float_t y        = slope*x + constant;
1036     Float_t theta    = TMath::ATan2(x,y);
1037
1038     if(x < fCorner1.X() && y > fCorner1.Y()) {
1039       arclength = radius*theta;                        //One sector since theta is by definition half-hybrid
1040
1041     }
1042
1043   }
1044
1045   return arclength;
1046
1047
1048 }
1049 //____________________________________________________________________
1050 Float_t
1051 AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)
1052 {
1053   //Get length of strip assuming that corners are not cut away
1054   Float_t rad             = GetMaxR(ring)-GetMinR(ring);
1055   Float_t nStrips         = (ring == 'I' ? 512 : 256);
1056   Float_t nSec            = (ring == 'I' ? 20 : 40);
1057   Float_t segment         = rad / nStrips;
1058   Float_t basearc         = 2*TMath::Pi() / (0.5*nSec); // One hybrid: 36 degrees inner, 18 outer
1059   Float_t radius          = GetMinR(ring) + strip*segment;
1060   Float_t basearclength   = 0.5*basearc * radius;                // One sector
1061
1062   return basearclength;
1063 }
1064 //____________________________________________________________________
1065 Int_t    AliFMDAnaParameters::GetFirstEtaBinFromMap(Int_t vtxbin, Int_t det, Char_t ring)
1066 {
1067   TH2F* hBg = GetBackgroundCorrection(det,ring,vtxbin);
1068
1069   if(det == 0) return hBg->GetXaxis()->FindBin(-1.95);
1070
1071   Int_t firstbin = -1;
1072   Int_t nNonZeroFirst = 0;
1073
1074   for(Int_t i=1;i<=hBg->GetNbinsX();i++) {
1075     if(nNonZeroFirst == fNumberOfEtaBinsToCut && firstbin==-1) firstbin = i;
1076
1077     for(Int_t j=1;j<=hBg->GetNbinsY();j++) {
1078
1079       Float_t value = hBg->GetBinContent(i,j);
1080
1081       if(value > 0.001 && nNonZeroFirst<fNumberOfEtaBinsToCut)
1082         {nNonZeroFirst++; break;}
1083
1084
1085     }
1086   }
1087
1088   return firstbin;
1089
1090 }
1091 //____________________________________________________________________
1092 Int_t    AliFMDAnaParameters::GetLastEtaBinFromMap(Int_t vtxbin, Int_t det, Char_t ring)
1093 {
1094   TH2F* hBg = GetBackgroundCorrection(det,ring,vtxbin);
1095   Int_t lastbin=-1;
1096   Int_t nNonZeroLast = 0;
1097
1098   if(det == 0) return hBg->GetXaxis()->FindBin(1.95);
1099
1100   for(Int_t i=hBg->GetNbinsX();i>0;i--) {
1101     if(nNonZeroLast == fNumberOfEtaBinsToCut && lastbin==-1) lastbin = i;
1102
1103     for(Int_t j=1;j<=hBg->GetNbinsY();j++) {
1104
1105       Float_t value = hBg->GetBinContent(i,j);
1106
1107       if(value > 0.001 && nNonZeroLast<fNumberOfEtaBinsToCut)
1108         {nNonZeroLast++; break; }
1109
1110
1111     }
1112   }
1113
1114   return lastbin;
1115 }
1116
1117 //____________________________________________________________________
1118 Int_t    AliFMDAnaParameters::GetFirstEtaBinToInclude(Int_t vtxbin, Int_t det, Char_t ring)
1119 {
1120   Int_t ringNumber = (ring == 'I' ? 0 : 1);
1121   return fEtaLowBinLimits.GetBinContent(det,ringNumber,vtxbin);
1122
1123 }
1124
1125 //____________________________________________________________________
1126 Int_t    AliFMDAnaParameters::GetLastEtaBinToInclude(Int_t vtxbin, Int_t det, Char_t ring)
1127 {
1128   Int_t ringNumber = (ring == 'I' ? 0 : 1);
1129   return fEtaHighBinLimits.GetBinContent(det,ringNumber,vtxbin);
1130
1131 }
1132 //____________________________________________________________________
1133 //
1134 // EOF
1135 //