FMD/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 "AliESDEvent.h"
  36 #include "AliESDVertex.h"
  37
  38 //====================================================================
  39 ClassImp(AliFMDAnaParameters)
  40 #if 0
  41   ; // This is here to keep Emacs for indenting the next line
  42 #endif
  43
  44 //const char* AliFMDAnaParameters::fgkBackgroundCorrection  = "FMD/Correction/Background";
  45 //const char* AliFMDAnaParameters::fgkEnergyDists = "FMD/Correction/EnergyDistribution";
  46 const char* AliFMDAnaParameters::fgkBackgroundID = "background";
  47 const char* AliFMDAnaParameters::fgkEnergyDistributionID = "energydistributions";
  48 const char* AliFMDAnaParameters::fgkEventSelectionEffID  = "eventselectionefficiency";
  49 //____________________________________________________________________
  50 AliFMDAnaParameters* AliFMDAnaParameters::fgInstance = 0;
  51
  52 //____________________________________________________________________
  53
  54 AliFMDAnaParameters*
  55 AliFMDAnaParameters::Instance()
  56 {
  57   // Get static instance
  58   if (!fgInstance) fgInstance = new AliFMDAnaParameters;
  59   return fgInstance;
  60 }
  61
  62 //____________________________________________________________________
  63 AliFMDAnaParameters::AliFMDAnaParameters() :
  64   fIsInit(kFALSE),
  65   fBackground(0),
  66   fEnergyDistribution(0),
  67   fEventSelectionEfficiency(0),
  68   fCorner1(4.2231, 26.6638),
  69   fCorner2(1.8357, 27.9500),
  70   fEnergyPath("$ALICE_ROOT/FMD/Correction/EnergyDistribution"),
  71   fBackgroundPath("$ALICE_ROOT/FMD/Correction/Background"),
  72   fEventSelectionEffPath("$ALICE_ROOT/FMD/Correction/EventSelectionEfficiency"),
  73   fProcessPrimary(kFALSE),
  74   fProcessHits(kFALSE),
  75   fTrigger(kMB1),
  76   fEnergy(k10000),
  77   fMagField(k5G),
  78   fSpecies(kPP)
  79 {
  80
  81
  82   //fVerticies.Add(new TVector2(4.2231, 26.6638));
  83   // fVerticies.Add(new TVector2(1.8357, 27.9500));
  84   // Default constructor
  85 }
  86 //____________________________________________________________________
  87 char* AliFMDAnaParameters::GetPath(const char* species) {
  88
  89   char* path ;
  90
  91   if(species == fgkBackgroundID)
  92     path = Form("%s/%s_%d_%d_%d_%d_%d_%d.root",
  93                 fBackgroundPath.Data(),
  94                 fgkBackgroundID,
  95                 fEnergy,
  96                 fTrigger,
  97                 fMagField,
  98                 fSpecies,
  99                 0,
 100                 0);
 101   if(species == fgkEnergyDistributionID)
 102     path = Form("%s/%s_%d_%d_%d_%d_%d_%d.root",
 103                 fEnergyPath.Data(),
 104                 fgkEnergyDistributionID,
 105                 fEnergy,
 106                 fTrigger,
 107                 fMagField,
 108                 fSpecies,
 109                 0,
 110                 0);
 111   if(species == fgkEventSelectionEffID)
 112     path = Form("%s/%s_%d_%d_%d_%d_%d_%d.root",
 113                 fEventSelectionEffPath.Data(),
 114                 fgkEventSelectionEffID,
 115                 fEnergy,
 116                 fTrigger,
 117                 fMagField,
 118                 fSpecies,
 119                 0,
 120                 0);
 121
 122   return path;
 123 }
 124 //____________________________________________________________________
 125 void AliFMDAnaParameters::Init(Bool_t forceReInit, UInt_t what)
 126 {
 127   // Initialize the parameters manager.  We need to get stuff from the
 128   // CDB here.
 129   if (forceReInit) fIsInit = kFALSE;
 130   if (fIsInit) return;
 131   if (what & kBackgroundCorrection)       InitBackground();
 132   if (what & kEnergyDistributions)        InitEnergyDists();
 133   if (what & kEventSelectionEfficiency)   InitEventSelectionEff();
 134
 135   fIsInit = kTRUE;
 136 }
 137 //____________________________________________________________________
 138
 139 void AliFMDAnaParameters::InitBackground() {
 140
 141   //AliCDBEntry*   background = GetEntry(fgkBackgroundCorrection);
 142
 143   TFile* fin = TFile::Open(GetPath(fgkBackgroundID));
 144
 145   if (!fin) return;
 146
 147   fBackground = dynamic_cast<AliFMDAnaCalibBackgroundCorrection*>(fin->Get(fgkBackgroundID));
 148   if (!fBackground) AliFatal("Invalid background object from CDB");
 149
 150 }
 151
 152 //____________________________________________________________________
 153
 154 void AliFMDAnaParameters::InitEnergyDists() {
 155
 156   TFile* fin = TFile::Open(GetPath(fgkEnergyDistributionID));
 157   //AliCDBEntry*   edist = GetEntry(fgkEnergyDists);
 158   if (!fin) return;
 159
 160   fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(fin->Get(fgkEnergyDistributionID));
 161
 162   if (!fEnergyDistribution) AliFatal("Invalid background object from CDB");
 163
 164 }
 165
 166 //____________________________________________________________________
 167
 168 void AliFMDAnaParameters::InitEventSelectionEff() {
 169
 170   //AliCDBEntry*   background = GetEntry(fgkBackgroundCorrection);
 171   TFile* fin = TFile::Open(GetPath(fgkEventSelectionEffID));
 172
 173   if (!fin) return;
 174
 175   fEventSelectionEfficiency = dynamic_cast<AliFMDAnaCalibEventSelectionEfficiency*>(fin->Get(fgkEventSelectionEffID));
 176   if (!fEventSelectionEfficiency) AliFatal("Invalid background object from CDB");
 177
 178 }
 179 //____________________________________________________________________
 180 Float_t AliFMDAnaParameters::GetVtxCutZ() {
 181
 182   if(!fIsInit) {
 183     AliWarning("Not initialized yet. Call Init() to remedy");
 184     return -1;
 185   }
 186
 187   return fBackground->GetVtxCutZ();
 188 }
 189
 190 //____________________________________________________________________
 191 Int_t AliFMDAnaParameters::GetNvtxBins() {
 192
 193   if(!fIsInit) {
 194     AliWarning("Not initialized yet. Call Init() to remedy");
 195     return -1;
 196   }
 197
 198   return fBackground->GetNvtxBins();
 199 }
 200 //____________________________________________________________________
 201 TH1F* AliFMDAnaParameters::GetEnergyDistribution(Int_t det, Char_t ring, Float_t eta) {
 202
 203   return fEnergyDistribution->GetEnergyDistribution(det, ring, eta);
 204 }
 205 //____________________________________________________________________
 206 Float_t AliFMDAnaParameters::GetSigma(Int_t det, Char_t ring, Float_t eta) {
 207
 208   if(!fIsInit) {
 209     AliWarning("Not initialized yet. Call Init() to remedy");
 210     return 0;
 211   }
 212
 213   TH1F* hEnergyDist       = GetEnergyDistribution(det,ring, eta);
 214   TF1*  fitFunc           = hEnergyDist->GetFunction("FMDfitFunc");
 215   if(!fitFunc) {
 216     AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
 217     return 1024;
 218   }
 219   Float_t sigma           = fitFunc->GetParameter(2);
 220   return sigma;
 221 }
 222
 223
 224 //____________________________________________________________________
 225 Float_t AliFMDAnaParameters::GetMPV(Int_t det, Char_t ring, Float_t eta) {
 226
 227   if(!fIsInit) {
 228     AliWarning("Not initialized yet. Call Init() to remedy");
 229     return 0;
 230   }
 231
 232   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 233   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 234   if(!fitFunc) {
 235     AliWarning(Form("No function for FMD%d%c, eta %f",det,ring,eta));
 236     return 1024;
 237   }
 238
 239   Float_t mpv           = fitFunc->GetParameter(1);
 240   return mpv;
 241 }
 242 //____________________________________________________________________
 243 Float_t AliFMDAnaParameters::Get2MIPWeight(Int_t det, Char_t ring, Float_t eta) {
 244
 245   if(!fIsInit) {
 246     AliWarning("Not initialized yet. Call Init() to remedy");
 247     return 0;
 248   }
 249
 250   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 251   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 252   if(!fitFunc) return 0;
 253   Float_t twoMIPweight    = fitFunc->GetParameter(3);
 254
 255
 256
 257   if(twoMIPweight < 1e-05)
 258     twoMIPweight = 0;
 259
 260   return twoMIPweight;
 261 }
 262 //____________________________________________________________________
 263 Float_t AliFMDAnaParameters::Get3MIPWeight(Int_t det, Char_t ring, Float_t eta) {
 264
 265   if(!fIsInit) {
 266     AliWarning("Not initialized yet. Call Init() to remedy");
 267     return 0;
 268   }
 269
 270   TH1F* hEnergyDist     = GetEnergyDistribution(det,ring,eta);
 271   TF1*  fitFunc         = hEnergyDist->GetFunction("FMDfitFunc");
 272   if(!fitFunc) return 0;
 273   Float_t threeMIPweight    = fitFunc->GetParameter(4);
 274
 275   if(threeMIPweight < 1e-05)
 276     threeMIPweight = 0;
 277
 278   Float_t twoMIPweight    = fitFunc->GetParameter(3);
 279
 280   if(twoMIPweight < 1e-05)
 281     threeMIPweight = 0;
 282
 283   return threeMIPweight;
 284 }
 285 //____________________________________________________________________
 286 Int_t AliFMDAnaParameters::GetNetaBins() {
 287   return GetBackgroundCorrection(1,'I',0)->GetNbinsX();
 288
 289 }
 290 //____________________________________________________________________
 291 Float_t AliFMDAnaParameters::GetEtaMin() {
 292
 293   return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmin();
 294 }
 295 //____________________________________________________________________
 296 Float_t AliFMDAnaParameters::GetEtaMax() {
 297
 298 return GetBackgroundCorrection(1,'I',0)->GetXaxis()->GetXmax();
 299
 300 }
 301 //____________________________________________________________________
 302
 303 TH2F* AliFMDAnaParameters::GetBackgroundCorrection(Int_t det,
 304                                                    Char_t ring,
 305                                                    Int_t vtxbin) {
 306
 307   if(!fIsInit) {
 308     AliWarning("Not initialized yet. Call Init() to remedy");
 309     return 0;
 310   }
 311
 312
 313
 314   if(vtxbin > fBackground->GetNvtxBins()) {
 315     AliWarning(Form("No background object for vertex bin %d", vtxbin));
 316     return 0;
 317   }
 318
 319   return fBackground->GetBgCorrection(det,ring,vtxbin);
 320 }
 321 //____________________________________________________________________
 322
 323 TH1F* AliFMDAnaParameters::GetDoubleHitCorrection(Int_t det,
 324                                                   Char_t ring) {
 325
 326   if(!fIsInit) {
 327     AliWarning("Not initialized yet. Call Init() to remedy");
 328     return 0;
 329   }
 330
 331   return fBackground->GetDoubleHitCorrection(det,ring);
 332 }
 333 //_____________________________________________________________________
 334 Float_t AliFMDAnaParameters::GetEventSelectionEfficiency(Int_t vtxbin) {
 335   if(!fIsInit) {
 336     AliWarning("Not initialized yet. Call Init() to remedy");
 337     return 0;
 338   }
 339   return fEventSelectionEfficiency->GetCorrection(vtxbin);
 340
 341 }
 342 //_____________________________________________________________________
 343 Float_t AliFMDAnaParameters::GetMaxR(Char_t ring) const{
 344   Float_t radius = 0;
 345   if(ring == 'I')
 346     radius = 17.2;
 347   else if(ring == 'O')
 348     radius = 28.0;
 349   else
 350     AliWarning("Unknown ring - must be I or O!");
 351
 352   return radius;
 353 }
 354 //_____________________________________________________________________
 355 Float_t AliFMDAnaParameters::GetMinR(Char_t ring) const{
 356   Float_t radius = 0;
 357   if(ring == 'I')
 358     radius = 4.5213;
 359   else if(ring == 'O')
 360     radius = 15.4;
 361   else
 362     AliWarning("Unknown ring - must be I or O!");
 363
 364   return radius;
 365
 366 }
 367 //_____________________________________________________________________
 368 void AliFMDAnaParameters::SetCorners(Char_t ring) {
 369
 370   if(ring == 'I') {
 371     fCorner1.Set(4.9895, 15.3560);
 372     fCorner2.Set(1.8007, 17.2000);
 373   }
 374   else {
 375     fCorner1.Set(4.2231, 26.6638);
 376     fCorner2.Set(1.8357, 27.9500);
 377   }
 378
 379 }
 380 //_____________________________________________________________________
 381 Float_t AliFMDAnaParameters::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec) const
 382 {
 383   Int_t nsec = (ring == 'I' ? 20 : 40);
 384   Float_t basephi = 0;
 385   if(det == 1)
 386     basephi = 1.72787594;
 387   if(det == 2 && ring == 'I')
 388     basephi = 0.15707963;
 389   if(det == 2 && ring == 'O')
 390     basephi = 0.078539818;
 391   if(det == 3 && ring == 'I')
 392     basephi = 2.984513044;
 393   if(det == 3 && ring == 'O')
 394     basephi = 3.06305289;
 395
 396   Float_t step = 2*TMath::Pi() / nsec;
 397   Float_t phi = 0;
 398   if(det == 3)
 399     phi = basephi - sec*step;
 400   else
 401     phi = basephi + sec*step;
 402
 403   if(phi < 0)
 404     phi = phi +2*TMath::Pi();
 405   if(phi > 2*TMath::Pi() )
 406     phi = phi - 2*TMath::Pi();
 407
 408   return phi;
 409 }
 410 //_____________________________________________________________________
 411 Float_t AliFMDAnaParameters::GetEtaFromStrip(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip, Float_t zvtx) const
 412 {
 413   // AliFMDRing fmdring(ring);
 414   // fmdring.Init();
 415   Float_t   rad       = GetMaxR(ring)-GetMinR(ring);
 416   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
 417   Float_t   segment   = rad / nStrips;
 418   Float_t   r         = GetMinR(ring) + segment*strip;
 419   Float_t   z         = 0;
 420   Int_t hybrid = sec / 2;
 421
 422   if(det == 1) {
 423     if(!(hybrid%2)) z = 320.266; else z = 319.766;
 424   }
 425   if(det == 2 && ring == 'I' ) {
 426     if(!(hybrid%2)) z = 83.666; else z = 83.166;
 427   }
 428   if(det == 2 && ring == 'O' ) {
 429     if(!(hybrid%2)) z = 74.966; else z = 75.466;
 430   }
 431   if(det == 3 && ring == 'I' ) {
 432     if(!(hybrid%2)) z = -63.066; else z = -62.566;
 433   }
 434   if(det == 3 && ring == 'O' ) {
 435     if(!(hybrid%2)) z = -74.966; else z = -75.466;
 436   }
 437
 438   //std::cout<<det<<"   "<<ring<<"   "<<sec<<"   "<<hybrid<<"    "<<z<<std::endl;
 439
 440   // Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
 441   Float_t   theta = TMath::ATan2(r,z-zvtx);
 442   Float_t   eta   = -1*TMath::Log(TMath::Tan(0.5*theta));
 443
 444   return eta;
 445 }
 446
 447 //_____________________________________________________________________
 448
 449 void AliFMDAnaParameters::GetVertex(AliESDEvent* esd, Double_t* vertexXYZ)
 450 {
 451   const AliESDVertex* vertex = 0;
 452   vertex = esd->GetPrimaryVertex();
 453   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 454     vertex = esd->GetPrimaryVertexSPD();
 455   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 456     vertex = esd->GetPrimaryVertexTPC();
 457   if(!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
 458     vertex = esd->GetVertex();
 459   if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
 460     vertex->GetXYZ(vertexXYZ);
 461     //std::cout<<vertex->GetName()<<"   "<< vertex->GetTitle() <<"   "<< vertex->GetZv()<<std::endl;
 462     return;
 463   }
 464   else if (esd->GetESDTZERO()) {
 465     vertexXYZ[0] = 0;
 466     vertexXYZ[1] = 0;
 467     vertexXYZ[2] = esd->GetT0zVertex();
 468
 469     return;
 470   }
 471
 472   return;
 473
 474 }
 475
 476 //____________________________________________________________________
 477 Bool_t AliFMDAnaParameters::IsEventTriggered(AliESDEvent *esd) {
 478   // check if the event was triggered
 479   ULong64_t triggerMask = esd->GetTriggerMask();
 480
 481   // definitions from p-p.cfg
 482   ULong64_t spdFO = (1 << 14);
 483   ULong64_t v0left = (1 << 11);
 484   ULong64_t v0right = (1 << 12);
 485
 486   switch (fTrigger) {
 487   case kMB1: {
 488     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 489       return kTRUE;
 490     break;
 491   }
 492   case kMB2: {
 493     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 494       return kTRUE;
 495     break;
 496   }
 497   case kSPDFASTOR: {
 498     if (triggerMask & spdFO)
 499       return kTRUE;
 500     break;
 501   }
 502   case kTEST_NOCTP: {
 503     return kTRUE;
 504     break;
 505   }
 506   }//switch
 507
 508   return kFALSE;
 509 }
 510
 511 //____________________________________________________________________
 512 Float_t
 513 AliFMDAnaParameters::GetStripLength(Char_t ring, UShort_t strip)
 514 {
 515   //AliFMDRing fmdring(ring);
 516   // fmdring.Init();
 517
 518   Float_t rad        = GetMaxR(ring)-GetMinR(ring);
 519   Float_t   nStrips   = (ring == 'I' ? 512 : 256);
 520   Float_t segment    = rad / nStrips;
 521
 522   //TVector2* corner1  = fmdring.GetVertex(2);
 523   // TVector2* corner2  = fmdring.GetVertex(3);
 524
 525   SetCorners(ring);
 526   /*
 527   std::cout<<GetMaxR(ring)<<"   "<<fmdring.GetMaxR()<<std::endl;
 528   std::cout<<GetMinR(ring)<<"   "<<fmdring.GetMinR()<<std::endl;
 529   std::cout<<corner1->X()<<"   "<<fCorner1.X()<<std::endl;
 530   std::cout<<corner2->X()<<"   "<<fCorner2.X()<<std::endl;
 531   std::cout<<corner1->Y()<<"   "<<fCorner1.Y()<<std::endl;
 532   std::cout<<corner2->Y()<<"   "<<fCorner2.Y()<<std::endl;*/
 533   Float_t slope      = (fCorner1.Y() - fCorner2.Y()) / (fCorner1.X() - fCorner2.X());
 534   Float_t constant   = (fCorner2.Y()*fCorner1.X()-(fCorner2.X()*fCorner1.Y())) / (fCorner1.X() - fCorner2.X());
 535   Float_t radius     = GetMinR(ring) + strip*segment;
 536
 537   Float_t d          = TMath::Power(TMath::Abs(radius*slope),2) + TMath::Power(radius,2) - TMath::Power(constant,2);
 538
 539   Float_t arclength  = GetBaseStripLength(ring,strip);
 540   if(d>0) {
 541
 542     Float_t x        = (-1*TMath::Sqrt(d) -slope*constant) / (1+TMath::Power(slope,2));
 543     Float_t y        = slope*x + constant;
 544     Float_t theta    = TMath::ATan2(x,y);
 545
 546     if(x < fCorner1.X() && y > fCorner1.Y()) {
 547       arclength = radius*theta;                        //One sector since theta is by definition half-hybrid
 548
 549     }
 550
 551   }
 552
 553   return arclength;
 554
 555
 556 }
 557 //____________________________________________________________________
 558 Float_t
 559 AliFMDAnaParameters::GetBaseStripLength(Char_t ring, UShort_t strip)
 560 {
 561   // AliFMDRing fmdring(ring);
 562   // fmdring.Init();
 563   Float_t rad             = GetMaxR(ring)-GetMinR(ring);
 564   Float_t nStrips         = (ring == 'I' ? 512 : 256);
 565   Float_t nSec            = (ring == 'I' ? 20 : 40);
 566   Float_t segment         = rad / nStrips;
 567   Float_t basearc         = 2*TMath::Pi() / (0.5*nSec); // One hybrid: 36 degrees inner, 18 outer
 568   Float_t radius          = GetMinR(ring) + strip*segment;
 569   Float_t basearclength   = 0.5*basearc * radius;                // One sector
 570
 571   return basearclength;
 572 }
 573 //____________________________________________________________________
 574 //
 575 // EOF
 576 //