PWG2/FORWARD/analysis/AliFMDAnaParameters.cxx

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