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