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
 380   if(!fitFunc) {
 381     AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
 382     std::cout<<hEnergyDist->GetEntries()<<"    "<<GetEtaBin(eta)<<std::endl;
 383     return 1024;
 384   }
 385
 386   Float_t mpv           = fitFunc->GetParameter(1);
 387   return mpv;
 388 }
 389 //____________________________________________________________________
 390 Float_t AliFMDAnaParameters::GetConstant(Int_t det, Char_t ring, Float_t eta) {
 391   //Get constant parameter of Landau fits
 392   if(!fIsInit) {
 393     AliWarning("Not initialized yet. Call Init() to remedy");
 394     return 0;
 395   }
 396
 397   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 398   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 399   if(!fitFunc) {
 400     AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
 401     return 0;
 402   }
 403
 404   Float_t mpv           = fitFunc->GetParameter(0);
 405   return mpv;
 406 }
 407 //____________________________________________________________________
 408 Float_t AliFMDAnaParameters::Get2MIPWeight(Int_t det, Char_t ring, Float_t eta) {
 409   //Get 2 MIP weights of convoluted Landau fits
 410   if(!fIsInit) {
 411     AliWarning("Not initialized yet. Call Init() to remedy");
 412     return 0;
 413   }
 414
 415   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 416   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 417   if(!fitFunc) return 0;
 418   Float_t twoMIPweight    = fitFunc->GetParameter(3);
 419
 420
 421
 422   if(twoMIPweight < 1e-05)
 423     twoMIPweight = 0;
 424
 425   return twoMIPweight;
 426 }
 427 //____________________________________________________________________
 428 Float_t AliFMDAnaParameters::Get3MIPWeight(Int_t det, Char_t ring, Float_t eta) {
 429   //Get 3 MIP weights of convoluted Landau fits
 430   if(!fIsInit) {
 431     AliWarning("Not initialized yet. Call Init() to remedy");
 432     return 0;
 433   }
 434
 435   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 436   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 437   if(!fitFunc) return 0;
 438   Float_t threeMIPweight    = fitFunc->GetParameter(4);
 439
 440   if(threeMIPweight < 1e-05)
 441     threeMIPweight = 0;
 442
 443   Float_t twoMIPweight    = fitFunc->GetParameter(3);
 444
 445   if(twoMIPweight < 1e-05)
 446     threeMIPweight = 0;
 447
 448   return threeMIPweight;
 449 }
 450 //____________________________________________________________________
 451 Int_t AliFMDAnaParameters::GetNetaBins() {
 452   return GetBackgroundCorrection(1,'I',5)->GetNbinsX();
 453
 454 }
 455 //____________________________________________________________________
 456 Float_t AliFMDAnaParameters::GetEtaMin() {
 457
 458   return GetBackgroundCorrection(1,'I',5)->GetXaxis()->GetXmin();
 459 }
 460 //____________________________________________________________________
 461 Float_t AliFMDAnaParameters::GetEtaMax() {
 462
 463 return GetBackgroundCorrection(1,'I',5)->GetXaxis()->GetXmax();
 464
 465 }
 466 //____________________________________________________________________
 467 Int_t AliFMDAnaParameters::GetEtaBin(Float_t eta) {
 468
 469   TAxis testaxis(GetNetaBins(),GetEtaMin(),GetEtaMax());
 470   Int_t binnumber = testaxis.FindBin(eta) ;
 471
 472   return binnumber;
 473
 474 }
 475 //____________________________________________________________________
 476
 477 TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det,
 478                                                    Char_t ring,
 479                                                    Int_t vtxbin) {
 480   //Get background correction histogram
 481   if(!fIsInit) {
 482     AliWarning("Not initialized yet. Call Init() to remedy");
 483     return 0;
 484   }
 485
 486
 487
 488   if(vtxbin > fBackground->GetNvtxBins()) {
 489     AliWarning(Form("No background object for vertex bin %d", vtxbin));
 490     return 0;
 491   }
 492
 493   return fBackground->GetBgCorrection(det,ring,vtxbin);
 494 }
 495 //____________________________________________________________________
 496
 497 TH1F* AliFMDAnaParameters::GetDoubleHitCorrection(Int_t det,
 498                                                   Char_t ring) {
 499   //Get correction for several hits in strips for p+p
 500   if(!fIsInit) {
 501     AliWarning("Not initialized yet. Call Init() to remedy");
 502     return 0;
 503   }
 504
 505   return fBackground->GetDoubleHitCorrection(det,ring);
 506 }
 507 //_____________________________________________________________________
 508 TH1F* AliFMDAnaParameters::GetSPDDeadCorrection(Int_t vtxbin) {
 509
 510   //Get correction for several hits in strips for p+p
 511   if(!fIsInit) {
 512     AliWarning("Not initialized yet. Call Init() to remedy");
 513     return 0;
 514   }
 515
 516   return fBackground->GetSPDDeadCorrection(vtxbin);
 517 }
 518 //_____________________________________________________________________
 519 Float_t AliFMDAnaParameters::GetEventSelectionEfficiency(Int_t vtxbin) {
 520   //Get event selection efficiency object
 521   if(!fIsInit) {
 522     AliWarning("Not initialized yet. Call Init() to remedy");
 523     return 0;
 524   }
 525   return fEventSelectionEfficiency->GetCorrection(vtxbin);
 526
 527 }
 528 //_____________________________________________________________________
 529 Float_t AliFMDAnaParameters::GetVtxSelectionEffFromMC() {
 530
 531    if(!fIsInit) {
 532     AliWarning("Not initialized yet. Call Init() to remedy");
 533     return 0;
 534    }
 535    return fEventSelectionEfficiency->GetVtxToTriggerRatio();
 536
 537
 538 }
 539 //_____________________________________________________________________
 540 TH2F* AliFMDAnaParameters::GetEventSelectionEfficiency(Int_t vtxbin, Char_t ring) {
 541   //Get event selection efficiency object
 542   if(!fIsInit) {
 543     AliWarning("Not initialized yet. Call Init() to remedy");
 544     return 0;
 545   }
 546   return fEventSelectionEfficiency->GetCorrection(vtxbin,ring);
 547
 548 }
 549 //_____________________________________________________________________
 550 TH1F* AliFMDAnaParameters::GetSharingEfficiency(Int_t det, Char_t ring, Int_t vtxbin) {
 551   //Get sharing efficiency object
 552   if(!fIsInit) {
 553     AliWarning("Not initialized yet. Call Init() to remedy");
 554     return 0;
 555   }
 556
 557   return fSharingEfficiency->GetSharingEff(det,ring,vtxbin);
 558
 559 }
 560 //_____________________________________________________________________
 561 TH1F* AliFMDAnaParameters::GetSharingEfficiencyTrVtx(Int_t det, Char_t ring, Int_t vtxbin) {
 562   //Get sharing efficiency object TrVtx
 563   if(!fIsInit) {
 564     AliWarning("Not initialized yet. Call Init() to remedy");
 565     return 0;
 566   }
 567
 568   return fSharingEfficiency->GetSharingEffTrVtx(det,ring,vtxbin);
 569
 570 }
 571 //_____________________________________________________________________
 572 Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) const {
 573   //Get max R of ring
 574   Float_t radius = 0;
 575   if(ring == 'I')
 576     radius = 17.2;
 577   else if(ring == 'O')
 578     radius = 28.0;
 579   else
 580     AliWarning("Unknown ring - must be I or O!");
 581
 582   return radius;
 583 }
 584 //_____________________________________________________________________
 585 Float_t AliFMDAnaParameters::GetMinR(Char_t ring) const{
 586   //Get min R of ring
 587   Float_t radius = 0;
 588   if(ring == 'I')
 589     radius = 4.5213;
 590   else if(ring == 'O')
 591     radius = 15.4;
 592   else
 593     AliWarning("Unknown ring - must be I or O!");
 594
 595   return radius;
 596
 597 }
 598 //_____________________________________________________________________
 599 void AliFMDAnaParameters::SetCorners(Char_t ring) {
 600   //Set corners (taken from nominal geometry)
 601   if(ring == 'I') {
 602     fCorner1.Set(4.9895, 15.3560);
 603     fCorner2.Set(1.8007, 17.2000);
 604   }
 605   else {
 606     fCorner1.Set(4.2231, 26.6638);
 607     fCorner2.Set(1.8357, 27.9500);
 608   }
 609
 610 }
 611 //_____________________________________________________________________
 612 Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) const
 613 {
 614   //Get phi from sector
 615   Int_t nsec = (ring == 'I' ? 20 : 40);
 616   Float_t basephi = 0;
 617   if(det == 1)
 618     basephi = 1.72787594;
 619   if(det == 2 && ring == 'I')
 620     basephi = 0.15707963;
 621   if(det == 2 && ring == 'O')
 622     basephi = 0.078539818;
 623   if(det == 3 && ring == 'I')
 624     basephi = 2.984513044;
 625   if(det == 3 && ring == 'O')
 626     basephi = 3.06305289;
 627
 628   Float_t step = 2*TMath::Pi() / nsec;
 629   Float_t phi = 0;
 630   if(det == 3)
 631     phi = basephi - sec*step;
 632   else
 633     phi = basephi + sec*step;
 634
 635   if(phi < 0)
 636     phi = phi +2*TMath::Pi();
 637   if(phi > 2*TMath::Pi() )
 638     phi = phi - 2*TMath::Pi();
 639
 640   return phi;
 641 }
 642 //_____________________________________________________________________
 643 Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) const
 644 {
 645   //Calculate eta from strip with vertex (redundant with AliESDFMD::Eta)
 646   Float_t   rad       = GetMaxR(ring)-GetMinR(ring);
 647   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
 648   Float_t   segment   = rad / nStrips;
 649   Float_t   r         = GetMinR(ring) + segment*strip;
 650   Float_t   z         = 0;
 651   Int_t hybrid = sec / 2;
 652
 653   if(det == 1) {
 654     if(!(hybrid%2)) z = 320.266; else z = 319.766;
 655   }
 656   if(det == 2 && ring == 'I' ) {
 657     if(!(hybrid%2)) z = 83.666; else z = 83.166;
 658   }
 659   if(det == 2 && ring == 'O' ) {
 660     if(!(hybrid%2)) z = 74.966; else z = 75.466;
 661   }
 662   if(det == 3 && ring == 'I' ) {
 663     if(!(hybrid%2)) z = -63.066; else z = -62.566;
 664   }
 665   if(det == 3 && ring == 'O' ) {
 666     if(!(hybrid%2)) z = -74.966; else z = -75.466;
 667   }
 668
 669   //std::cout<<det<<"   "<<ring<<"   "<<sec<<"   "<<hybrid<<"    "<<z<<std::endl;
 670
 671   // Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 672   Float_t   theta = TMath::ATan2(r,z-zvtx);
 673   Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 674
 675   return eta;
 676 }
 677
 678 //_____________________________________________________________________
 679
 680 Bool_t AliFMDAnaParameters::GetVertex(const AliESDEvent* esd, Double_t* vertexXYZ)
 681 {
 682   //Get the vertex from the ESD
 683   const AliESDVertex* vertex = 0;
 684   vertex = esd->GetPrimaryVertexSPD();
 685
 686   if(vertex)
 687     vertex->GetXYZ(vertexXYZ);
 688
 689   //if(vertexXYZ[0] == 0 || vertexXYZ[1] == 0 )
 690   //  return kFALSE;
 691
 692   if(vertex->GetNContributors() <= 0)
 693     return kFALSE;
 694
 695   if(vertex->GetZRes() > 0.1 )
 696     return kFALSE;
 697
 698   return vertex->GetStatus();
 699
 700 }
 701 //____________________________________________________________________
 702 Bool_t AliFMDAnaParameters::IsEventTriggered(const AliESDEvent *esd, Trigger trig) {
 703   //Did we have trig trigger ?
 704   Trigger old = fTrigger;
 705   fTrigger = trig;
 706   Bool_t retval = IsEventTriggered(esd);
 707   fTrigger = old;
 708   return retval;
 709
 710 }
 711 //____________________________________________________________________
 712 Bool_t AliFMDAnaParameters::IsEventTriggered(const AliESDEvent *esd) const {
 713   // check if the event was triggered
 714
 715   if (fCentralSelection) return kTRUE;
 716   ULong64_t triggerMask = esd->GetTriggerMask();
 717
 718   TString triggers = esd->GetFiredTriggerClasses();
 719
 720   //if(triggers.Contains("CINT1B-ABCE-NOPF-ALL") || triggers.Contains("CINT1B-E-NOPF-ALL")) return kTRUE;
 721   //else return kFALSE;
 722   //if(triggers.Contains("CINT1B-E-NOPF-ALL"))    return kFALSE;
 723
 724   // if(triggers.Contains("CINT1A-ABCE-NOPF-ALL")) return kFALSE;
 725   // if(triggers.Contains("CINT1C-ABCE-NOPF-ALL")) return kFALSE;
 726
 727   // definitions from p-p.cfg
 728   ULong64_t spdFO = (1 << 14);
 729   ULong64_t v0left = (1 << 11);
 730   ULong64_t v0right = (1 << 12);
 731   AliTriggerAnalysis tAna;
 732
 733   //REMOVE WHEN FINISHED PLAYING WITH TRIGGERS!
 734   //fPhysicsSelection->IsCollisionCandidate(esd);
 735   if(!fRealData) {
 736     fPhysicsSelection->SetAnalyzeMC(kTRUE);
 737   }
 738   switch (fTrigger) {
 739   case kMB1: {
 740     // if(fRealData) {
 741       if( fPhysicsSelection->IsCollisionCandidate(esd))
 742         return kTRUE;
 743       //}
 744       //else {
 745       // if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 746       //        return kTRUE;
 747       break;
 748       //}
 749   }
 750   case kMB2: {
 751     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 752       return kTRUE;
 753     break;
 754   }
 755   case kSPDFASTOR: {
 756     if (triggerMask & spdFO)
 757       return kTRUE;
 758     break;
 759   }
 760   case kNOCTP: {
 761     return kTRUE;
 762     break;
 763   }
 764   case kEMPTY: {
 765     /*
 766     const AliMultiplicity* testmult = esd->GetMultiplicity();
 767     Int_t nTrackLets = testmult->GetNumberOfTracklets();
 768     Int_t nClusters  = testmult->GetNumberOfSingleClusters();
 769     Int_t fastOR = tAna.SPDFiredChips(esd, 0);
 770
 771     const AliESDVertex* vertex = 0;
 772     vertex = esd->GetPrimaryVertexSPD();
 773     Bool_t vtxstatus = vertex->GetStatus();
 774     if(vertex->GetNContributors() <= 0)
 775       vtxstatus = kFALSE;
 776
 777     if(vertex->GetZRes() > 0.1 )
 778       vtxstatus = kFALSE;
 779
 780     Bool_t v0ABG = tAna.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0ABG);
 781
 782     Bool_t v0CBG = tAna.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0CBG);
 783     Bool_t v0A = tAna.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0A);
 784     Bool_t v0C = tAna.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kV0C);*/
 785     //    if(triggers.Contains("CINT1A-ABCE-NOPF-ALL") || triggers.Contains("CINT1C-ABCE-NOPF-ALL"))
 786     if(triggers.Contains("CBEAMB-ABCE-NOPF-ALL"))
 787       return kTRUE;
 788     break;
 789   }
 790
 791   }//switch
 792
 793   return kFALSE;
 794 }
 795
 796 //____________________________________________________________________
 797 Float_t
 798 AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)
 799 {
 800   //Get length of a strip
 801
 802   Float_t rad        = GetMaxR(ring)-GetMinR(ring);
 803   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
 804   Float_t segment    = rad / nStrips;
 805
 806   //TVector2* corner1  = fmdring.GetVertex(2);
 807   // TVector2* corner2  = fmdring.GetVertex(3);
 808
 809   SetCorners(ring);
 810   /*
 811   std::cout<<GetMaxR(ring)<<"   "<<fmdring.GetMaxR()<<std::endl;
 812   std::cout<<GetMinR(ring)<<"   "<<fmdring.GetMinR()<<std::endl;
 813   std::cout<<corner1->X()<<"   "<<fCorner1.X()<<std::endl;
 814   std::cout<<corner2->X()<<"   "<<fCorner2.X()<<std::endl;
 815   std::cout<<corner1->Y()<<"   "<<fCorner1.Y()<<std::endl;
 816   std::cout<<corner2->Y()<<"   "<<fCorner2.Y()<<std::endl;*/
 817   Float_t slope      = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
 818   Float_t constant   = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
 819   Float_t radius     = GetMinR(ring) + strip*segment;
 820
 821   Float_t d          = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
 822
 823   Float_t arclength  = GetBaseStripLength(ring,strip);
 824   if(d>0) {
 825
 826     Float_t x        = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
 827     Float_t y        = slope*x + constant;
 828     Float_t theta    = TMath::ATan2(x,y);
 829
 830     if(x < fCorner1.X() && y > fCorner1.Y()) {
 831       arclength = radius*theta;                        //One sector since theta is by definition half-hybrid
 832
 833     }
 834
 835   }
 836
 837   return arclength;
 838
 839
 840 }
 841 //____________________________________________________________________
 842 Float_t
 843 AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)
 844 {
 845   //Get length of strip assuming that corners are not cut away
 846   Float_t rad             = GetMaxR(ring)-GetMinR(ring);
 847   Float_t nStrips         = (ring == 'I' ? 512 : 256);
 848   Float_t nSec            = (ring == 'I' ? 20 : 40);
 849   Float_t segment         = rad / nStrips;
 850   Float_t basearc         = 2*TMath::Pi() / (0.5*nSec); // One hybrid: 36 degrees inner, 18 outer
 851   Float_t radius          = GetMinR(ring) + strip*segment;
 852   Float_t basearclength   = 0.5*basearc * radius;                // One sector
 853
 854   return basearclength;
 855 }
 856 //____________________________________________________________________
 857 Int_t    AliFMDAnaParameters::GetFirstEtaBinFromMap(Int_t vtxbin, Int_t det, Char_t ring)
 858 {
 859   TH2F* hBg = GetBackgroundCorrection(det,ring,vtxbin);
 860
 861   if(det == 0) return hBg->GetXaxis()->FindBin(-1.95);
 862
 863   Int_t firstbin = -1;
 864   Int_t nNonZeroFirst = 0;
 865
 866   for(Int_t i=1;i<=hBg->GetNbinsX();i++) {
 867     if(nNonZeroFirst == fNumberOfEtaBinsToCut && firstbin==-1) firstbin = i;
 868
 869     for(Int_t j=1;j<=hBg->GetNbinsY();j++) {
 870
 871       Float_t value = hBg->GetBinContent(i,j);
 872
 873       if(value > 0.001 && nNonZeroFirst<fNumberOfEtaBinsToCut)
 874         {nNonZeroFirst++; break;}
 875
 876
 877     }
 878   }
 879
 880   return firstbin;
 881
 882 }
 883 //____________________________________________________________________
 884 Int_t    AliFMDAnaParameters::GetLastEtaBinFromMap(Int_t vtxbin, Int_t det, Char_t ring)
 885 {
 886   TH2F* hBg = GetBackgroundCorrection(det,ring,vtxbin);
 887   Int_t lastbin=-1;
 888   Int_t nNonZeroLast = 0;
 889
 890   if(det == 0) return hBg->GetXaxis()->FindBin(1.95);
 891
 892   for(Int_t i=hBg->GetNbinsX();i>0;i--) {
 893     if(nNonZeroLast == fNumberOfEtaBinsToCut && lastbin==-1) lastbin = i;
 894
 895     for(Int_t j=1;j<=hBg->GetNbinsY();j++) {
 896
 897       Float_t value = hBg->GetBinContent(i,j);
 898
 899       if(value > 0.001 && nNonZeroLast<fNumberOfEtaBinsToCut)
 900         {nNonZeroLast++; break; }
 901
 902
 903     }
 904   }
 905
 906   return lastbin;
 907 }
 908
 909 //____________________________________________________________________
 910 Int_t    AliFMDAnaParameters::GetFirstEtaBinToInclude(Int_t vtxbin, Int_t det, Char_t ring)
 911 {
 912   Int_t ringNumber = (ring == 'I' ? 0 : 1);
 913   return fEtaLowBinLimits.GetBinContent(det,ringNumber,vtxbin);
 914
 915 }
 916
 917 //____________________________________________________________________
 918 Int_t    AliFMDAnaParameters::GetLastEtaBinToInclude(Int_t vtxbin, Int_t det, Char_t ring)
 919 {
 920   Int_t ringNumber = (ring == 'I' ? 0 : 1);
 921   return fEtaHighBinLimits.GetBinContent(det,ringNumber,vtxbin);
 922
 923 }
 924 //____________________________________________________________________
 925 //
 926 // EOF
 927 //