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