PWGLF/FORWARD/analysis2/AliDisplacedVertexSelection.cxx

   1 #include "AliDisplacedVertexSelection.h"
   2 // #include "AliAnalysisManager.h"
   3 #include "AliMCEvent.h"
   4 #include "AliMultiplicity.h"
   5 // #include "AliMCEventHandler.h"
   6 #include "AliHeader.h"
   7 #include "AliGenEventHeader.h"
   8 #include <iostream>
   9 #include <TROOT.h>
  10 #include <TH1D.h>
  11 #include <TH2D.h>
  12 #include "AliESDEvent.h"
  13 #include "AliESDZDC.h"
  14 ClassImp(AliDisplacedVertexSelection)
  15 #if 0
  16 ; // This is for Emacs
  17 #endif
  18
  19 //____________________________________________________________________
  20 AliDisplacedVertexSelection::AliDisplacedVertexSelection()
  21   : TObject(),
  22     fDeltaTdc(0),
  23     fSumTdc(0),
  24     fZdcEnergy(0),
  25     fZemEnergy(0),
  26     fCorrelationZemZdc(0),
  27     fCorrelationSumDelta(0),
  28     fVertexZ(kInvalidVtxZ),
  29     fCent(100),
  30     fHVertexZ(0),
  31     fHCent(0),
  32     fMC(kFALSE)
  33 {
  34 }
  35 //____________________________________________________________________
  36 AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
  37   : TObject(o),
  38     fDeltaTdc(o.fDeltaTdc),
  39     fSumTdc(o.fSumTdc),
  40     fZdcEnergy(o.fZdcEnergy),
  41     fZemEnergy(o.fZemEnergy),
  42     fCorrelationZemZdc(o.fCorrelationZemZdc),
  43     fCorrelationSumDelta(o.fCorrelationSumDelta),
  44     fVertexZ(kInvalidVtxZ),
  45     fCent(100),
  46     fHVertexZ(0),
  47     fHCent(0),
  48     fMC(kFALSE)
  49 {
  50 }
  51 //____________________________________________________________________
  52 AliDisplacedVertexSelection&
  53 AliDisplacedVertexSelection::operator=(const AliDisplacedVertexSelection& o)
  54 {
  55   if (&o == this) return *this;
  56
  57   fDeltaTdc  = o.fDeltaTdc;
  58   fSumTdc    = o.fSumTdc;
  59   fZdcEnergy = o.fZdcEnergy;
  60   fZemEnergy = o.fZemEnergy;
  61   fCorrelationZemZdc = o.fCorrelationZemZdc;
  62   fCorrelationSumDelta = o.fCorrelationSumDelta;
  63   fMC = o.fMC;
  64
  65   return *this;
  66 }
  67
  68 //____________________________________________________________________
  69 void
  70 AliDisplacedVertexSelection::SetupForData(TList* l,
  71                                           const char* /* name*/,
  72                                           Bool_t mc)
  73 {
  74   fMC = mc;
  75
  76   TList* out = new TList;
  77   out->SetName("displacedVertex");
  78   out->SetOwner();
  79   l->Add(out);
  80
  81   Double_t dVz   = 37.5;
  82   Double_t vzMin = (-kMaxK-.5) * dVz;
  83   Double_t vzMax = (+kMaxK+.5) * dVz;
  84
  85   fHVertexZ = new TH1D("vertexZ", "Interaction point Z",
  86                        2*kMaxK+1, vzMin, vzMax);
  87   fHVertexZ->SetXTitle("IP_{z} [cm]");
  88   fHVertexZ->SetYTitle("events");
  89   fHVertexZ->SetDirectory(0);
  90   fHVertexZ->SetFillColor(kRed+1);
  91   fHVertexZ->SetFillStyle(3001);
  92   out->Add(fHVertexZ);
  93
  94   Int_t    nCent   = 6;
  95   Double_t bCent[] = { 0, 5, 10, 20, 30, 40, 100 };
  96   fHCent = new TH1D("cent", "Centrality", nCent, bCent);
  97   fHCent->SetXTitle("Centrality [%]");
  98   fHCent->SetYTitle("events");
  99   fHCent->SetDirectory(0);
 100   fHCent->SetFillColor(kBlue+1);
 101   fHCent->SetFillStyle(3001);
 102   out->Add(fHCent);
 103
 104   Int_t    nDeltaT   = 2000;
 105   Double_t minDeltaT = -250;
 106   Double_t maxDeltaT = +250;
 107   fDeltaTdc          = new TH1D("DeltaTdc","#DeltaTDC",
 108                                 nDeltaT,minDeltaT,maxDeltaT);
 109   fDeltaTdc->SetXTitle("#DeltaTDC");
 110   fDeltaTdc->SetDirectory(0);
 111   out->Add(fDeltaTdc);
 112
 113   Int_t    nSumT   = nDeltaT;
 114   Double_t minSumT = minDeltaT;
 115   Double_t maxSumT = maxDeltaT;
 116   fSumTdc          = new TH1D("SumTdc","#sumTDC",nSumT,minSumT,maxSumT);
 117   fSumTdc->SetXTitle("#sumTDC");
 118   fSumTdc->SetDirectory(0);
 119   out->Add(fSumTdc);
 120
 121   Int_t    nZdcE   = 10000;
 122   Double_t minZdcE = 0;
 123   Double_t maxZdcE = 200000;
 124   fZdcEnergy           = new TH1D("ZDCEnergy","E_{ZDC}",nZdcE,minZdcE,maxZdcE);
 125   fZdcEnergy->SetXTitle("E_{ZDC}");
 126   fZdcEnergy->SetDirectory(0);
 127   out->Add(fZdcEnergy);
 128
 129   Int_t    nZemE   = 1000;
 130   Double_t minZemE = -100;
 131   Double_t maxZemE = 2900;
 132   fZemEnergy           = new TH1D("ZEMEnergy","E_{ZEM}",nZemE,minZemE,maxZemE);
 133   fZemEnergy->SetXTitle("E_{ZEM}");
 134   fZemEnergy->SetDirectory(0);
 135   out->Add(fZemEnergy);
 136
 137   fCorrelationZemZdc   = new TH2D("corrZEMZDC","E_{ZEM} vs E_{ZDC}",
 138                                   nZemE, minZemE, maxZemE,
 139                                   nZdcE, minZdcE, maxZdcE);
 140   fCorrelationZemZdc->SetXTitle("E_{ZEM}");
 141   fCorrelationZemZdc->SetYTitle("E_{ZDC}");
 142   fCorrelationZemZdc->SetDirectory(0);
 143   out->Add(fCorrelationZemZdc);
 144
 145   fCorrelationSumDelta =  new TH2D("corrSumDelta","#sumTDC vs #DeltaTDC",
 146                                    nSumT, minSumT, maxSumT,
 147                                    nDeltaT, minDeltaT, maxDeltaT);
 148   fCorrelationSumDelta->SetXTitle("#sum TDC");
 149   fCorrelationSumDelta->SetYTitle("#DeltaTDC");
 150   fCorrelationSumDelta->SetDirectory(0);
 151   out->Add(fCorrelationSumDelta);
 152
 153 }
 154
 155
 156 //____________________________________________________________________
 157 void
 158 AliDisplacedVertexSelection::Print(Option_t*) const
 159 {
 160 #if 0
 161   char ind[gROOT->GetDirLevel()+1];
 162   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 163   ind[gROOT->GetDirLevel()] = '\0';
 164   std::cout << std::boolalpha
 165             << std::noboolalpha << std::endl;
 166 #endif
 167 }
 168
 169 //____________________________________________________________________
 170 Float_t
 171 AliDisplacedVertexSelection::GetZemCorr(Int_t k, Bool_t minusminus) const
 172 {
 173   if (-kMaxK > k || k > kMaxK) return 0;
 174
 175   // Corrections for magnetic fields
 176   const Float_t kPlusPlus[21]   = { 0.6840,0.7879,0.8722,
 177                                     0.9370,0.9837,1.0137,
 178                                     1.0292,1.0327,1.0271,
 179                                     1.0152,1.0000,0.9844,
 180                                     0.9714,0.9634,0.9626,
 181                                     0.9708,0.9891,1.0181,
 182                                     1.0574,1.1060,1.1617};
 183   const Float_t kMoinsMoins[21] = { 0.7082,0.8012,0.8809,
 184                                     0.9447,0.9916,1.0220,
 185                                     1.0372,1.0395,1.0318,
 186                                     1.0174,1.0000,0.9830,
 187                                     0.9699,0.9635,0.9662,
 188                                     0.9794,1.0034,1.0371,
 189                                     1.0782,1.1224,1.1634};
 190
 191   // MC specific code by Hans, is it used - why are ++ and -- the
 192   // same?
 193   const Float_t kPlusPlusMC[21]   = { 0.68400, 0.78790, 0.87220,
 194                                       0.93700, 0.98370, 1.08982,
 195                                       1.03518, 1.01534, 1.01840,
 196                                       1.01493, 1.00000, 0.970083,
 197                                       0.979705,0.960576,0.925394,
 198                                       0.92167, 0.971241,1.08338,
 199                                       1.23517, 1.39308, 1.53943 };
 200   const Float_t kMoinsMoinsMC[21] = { 0.68400, 0.78790, 0.87220,
 201                                       0.9370,  0.98370, 1.08982,
 202                                       1.03518, 1.01534, 1.0184,
 203                                       1.01493, 1.00000, 0.970083,
 204                                       0.979705,0.960576,0.925394,
 205                                       0.92167, 0.971241,1.08338,
 206                                       1.23517, 1.39308, 1.53943 };
 207   Int_t i = k+kMaxK;
 208   if (fMC)
 209     return minusminus ? kMoinsMoinsMC[i] : kPlusPlusMC[i];
 210   return minusminus ? kMoinsMoins[i] : kPlusPlus[i];
 211 }
 212
 213
 214 //____________________________________________________________________
 215 Bool_t
 216 AliDisplacedVertexSelection::CheckOutlier(Int_t ivtx,
 217                                           const AliESDEvent* esd) const
 218 {
 219   if (fMC) return false;
 220   if (ivtx <= 4) return false; // Not outliers
 221
 222   // --- Find outliers -----------------------------------------------
 223   AliESDVZERO* esdV0 = esd->GetVZEROData();
 224   Double_t nClusters = esd->GetMultiplicity()->GetNumberOfITSClusters(1);
 225
 226   // parameter from the fit of VZERO multiplicity vs N SPD cluster
 227   // distribution
 228   const Double_t slp[8][21] = {{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//0
 229                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 230                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000},
 231                                {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//1
 232                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 233                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000},
 234                                {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//2
 235                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 236                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000},
 237                                {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//3
 238                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 239                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000},
 240                                {0.000,0.000,0.000,0.000,0.000,0.155,0.173,//4
 241                                 0.238,0.348,0.384,0.209,0.362,0.433,0.435,
 242                                 0.421,0.400,0.403,0.398,0.407,0.394,0.369},
 243                                {0.000,0.000,0.000,0.000,0.000,0.356,0.402,//5
 244                                 0.504,0.650,0.643,0.331,0.543,0.640,0.642,
 245                                 0.614,0.591,0.587,0.564,0.572,0.524,0.576},
 246                                {0.000,0.000,0.000,0.000,0.000,0.386,0.500,//6
 247                                 0.692,0.856,0.810,0.390,0.619,0.713,0.708,
 248                                 0.670,0.633,0.622,0.593,0.587,0.538,0.587},
 249                                {0.000,0.000,0.000,0.000,0.000,0.353,0.454,//7
 250                                 0.697,0.944,0.941,0.444,0.670,0.746,0.736,
 251                                 0.683,0.642,0.619,0.583,0.576,0.535,0.581}};
 252   const Double_t cst[8][21] = {{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//0
 253                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 254                                 0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
 255                                {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//1
 256                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 257                                 0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
 258                                {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//2
 259                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 260                                 0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
 261                                {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//3
 262                                 0.000,0.000,0.000,0.000,0.000,0.000,0.000,
 263                                 0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
 264                                {0.000,0.000,0.000,0.000,0.000,51.67,54.79,//4
 265                                 53.51,23.49,24.60,28.47,27.35,24.47,23.93,
 266                                 13.162,12.91,5.305,4.007,-14.34,-9.16,-9.43},
 267                                {0.000,0.000,0.000,0.000,0.000,136.34,88.84,
 268                                 84.85,24.31,27.92,35.64,33.15,32.25,23.07,
 269                                 13.746,-6.05,-11.95,1.995,-24.91,-22.55,-37.86},
 270                                {0.000,0.000,0.000,0.000,0.000,380.42,204.22,
 271                                 123.89,34.60,29.26,32.49,30.94,19.67,7.760,
 272                                 1.1010,-6.01,-15.66,-3.16,-18.34,-17.33,-22.17},
 273                                {0.000,0.000,0.000,0.000,0.000,128.00,105.04,
 274                                 114.05,15.96,21.67,30.11,30.69,27.66,0.363,-0.511,
 275                                 -17.67,-22.40,-11.84,-25.65,-24.29,-24.46}};
 276
 277   Bool_t outlier = kFALSE;
 278   // loop over VZERO rings
 279   for (int iring = 4; iring < 8; iring++) {
 280     Double_t multRing = 0;
 281     for (int iCh=0; iCh<8; iCh++) {
 282       Int_t idx = iCh+iring*8;
 283       multRing += esdV0->GetMultiplicity(idx);
 284     }
 285
 286     // Tigher cut for z=0.  Looser cut for suffician statistics and
 287     // see saturation effects if present.
 288     Double_t upFac    = (ivtx == 10 ? .35 : .42);
 289     Double_t downFac  = (ivtx == 10 ? .20 : ivtx < 8 ? .42 : .26);
 290     Double_t slpUP    = slp[iring][ivtx]*(1.+upFac); //upper cut
 291     Double_t slpDOWN  = slp[iring][ivtx]*(1.-downFac); //lower cut
 292     Double_t constant = cst[iring][ivtx];
 293     //Cut is applied for rings and vertex of interest
 294     if ((slpDOWN * nClusters + constant) > multRing ||
 295         (slpUP   * nClusters + constant) < multRing ) {
 296         outlier = kTRUE;
 297     }
 298   }
 299
 300   return outlier;
 301   //--- End of outlier code -----------------------------------------
 302 }
 303 //____________________________________________________________________
 304 Bool_t
 305 AliDisplacedVertexSelection::ProcessMC(const AliMCEvent* mcevent)
 306 {
 307   if (!fMC) return true;
 308   if (!mcevent) return false;
 309
 310   AliMCEvent*        event     = const_cast<AliMCEvent*>(mcevent);
 311   AliHeader*         header    = event->Header();
 312   AliGenEventHeader* genHeader = header->GenEventHeader();
 313   TArrayF            vertex;
 314   genHeader->PrimaryVertex(vertex);
 315   Double_t zvtx           = vertex.At(2);
 316   Double_t intTime        = genHeader->InteractionTime();
 317   const Double_t kVtx[16] = {-187.5, -150., -112.5, -75., -37.5, 0.,
 318                              37.5,   75.,  112.5, 150., 187.5, 225.,
 319                              262.5,  300.,  337.5, 375.};
 320   Int_t vtxClass          = -1;
 321   Float_t SigmaVtx        = 5.3;
 322   Float_t SigmaTime       = 0.18;
 323
 324   for(Int_t iVtx=0; iVtx<16; ++iVtx) {
 325     // Do not do something for nominal
 326     if (iVtx == 5) continue;
 327     if ((zvtx >= kVtx[iVtx] - 3.54 * SigmaVtx &&
 328          zvtx <= kVtx[iVtx] + 3.54 * SigmaVtx) &&
 329         (intTime >= (iVtx * 1.25e-9 - 6.25e-9 - 4 * SigmaTime) &&
 330          intTime <= (iVtx * 1.25e-9 - 6.25e-9 + 4 * SigmaTime)))  {
 331       vtxClass = iVtx;
 332       // break here to not do more
 333       break;
 334     }
 335   }
 336   if(vtxClass > -1) {
 337     fVertexZ = kVtx[vtxClass];
 338     return true;
 339   }
 340   return false;
 341 }
 342 //____________________________________________________________________
 343 Bool_t
 344 AliDisplacedVertexSelection::Process(const AliESDEvent* esd)
 345 {
 346   fVertexZ = kInvalidVtxZ; // Default vertex value
 347   fCent    = 100;  // Default centrality value
 348
 349   // Some constants
 350   const Float_t kZDCrefSum        = -66.5;
 351   const Float_t kZDCrefDelta      =  -2.10;
 352   const Float_t kZDCsigmaSum      =   3.25;
 353   const Float_t kZDCsigmaDelta    =   2.25;
 354   const Float_t kZDCsigmaSumSat   =   2.50;
 355   const Float_t kZDCsigmaDeltaSat =   1.35;
 356
 357   // --- Get the ESD object ------------------------------------------
 358   AliESDZDC* esdZDC = esd->GetESDZDC();
 359   if (!esdZDC) {
 360     AliWarning("No ZDC ESD object!");
 361     return false;
 362   }
 363   Double_t currentL3   = esd->GetCurrentL3();
 364   Double_t currentDipo = esd->GetCurrentDip();
 365
 366   // --- ZDC and ZEM energy signal -----------------------------------
 367   Double_t zdcEn      = (esdZDC->GetZDCN1Energy()+
 368                          esdZDC->GetZDCP1Energy()+
 369                          esdZDC->GetZDCN2Energy()+
 370                          esdZDC->GetZDCP2Energy());
 371   Double_t zemEn      = (esdZDC->GetZDCEMEnergy(0)+
 372                          esdZDC->GetZDCEMEnergy(1))/8.;
 373
 374   // HHD/Maxime inclusion - MC check!
 375   if (fMC) {
 376     zdcEn      = (2.9 * esdZDC->GetZDCN1Energy() +
 377                   7.2 * esdZDC->GetZDCP1Energy() +
 378                   3   * esdZDC->GetZDCN2Energy() +
 379                   8.7 * esdZDC->GetZDCP2Energy());
 380     zemEn      = 0.57 * (esdZDC->GetZDCEMEnergy(0) +
 381                          esdZDC->GetZDCEMEnergy(1));
 382   }
 383
 384   // --- Calculate DeltaT and sumT -----------------------------------
 385   Double_t deltaTdc = 0;
 386   Double_t sumTdc   = 0;
 387
 388   for(Int_t i = 0; i < 4; ++i) {
 389     if(esdZDC->GetZDCTDCData(10,i) != 0) {
 390       Double_t  tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
 391                                     esdZDC->GetZDCTDCData(14,i));
 392       Double_t  tdcC       = esdZDC->GetZDCTDCCorrected(10,i);
 393       for(Int_t j = 0; j < 4; ++j) {
 394         if(esdZDC->GetZDCTDCData(12,j) != 0) {
 395           Double_t   tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
 396                                          esdZDC->GetZDCTDCData(14,j));
 397           Double_t   tdcA       = esdZDC->GetZDCTDCCorrected(12,j);
 398           if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
 399             deltaTdc = tdcC-tdcA;
 400             sumTdc   = tdcC+tdcA;
 401           }
 402           else {
 403             deltaTdc = tdcCnoCorr-tdcAnoCorr;
 404             sumTdc   = tdcCnoCorr+tdcAnoCorr;
 405           }
 406         }
 407       }
 408     }
 409   }
 410   fDeltaTdc->Fill(deltaTdc);
 411   fSumTdc->Fill(sumTdc);
 412   fCorrelationSumDelta->Fill(sumTdc, deltaTdc);
 413
 414   // --- Find the vertex ---------------------------------------------
 415   Int_t ivtx = -1;
 416   if(deltaTdc!=0. || sumTdc!=0.) {
 417     Double_t fillVz = kInvalidVtxZ;
 418     for (Int_t k = -kMaxK; k <= kMaxK; ++k) {
 419       Float_t zsat  = 2.55F * k;
 420       Float_t delta = (k == 0 ? kZDCsigmaDelta : kZDCsigmaDeltaSat);
 421       Float_t sum   = (k == 0 ? kZDCsigmaSum   : kZDCsigmaSumSat);
 422       Float_t dT    = deltaTdc - kZDCrefDelta - zsat;
 423       Float_t sT    = sumTdc  - kZDCrefSum - zsat;
 424       Float_t check = dT * dT / delta / delta + sT * sT / sum  / sum;
 425       if (check > 1.0) continue;
 426       if (k == 0) {
 427         fillVz = 0;
 428         continue;
 429       }
 430       // Set the vertex
 431       fVertexZ = 37.5 * k;
 432       fillVz   = fVertexZ;
 433       ivtx     = k + 10; // Used for outlier calculations
 434       // Correct zem energy
 435       if(currentDipo>0 && currentL3>0) zemEn /= GetZemCorr(k,false);
 436       if(currentDipo<0 && currentL3<0) zemEn /= GetZemCorr(k,true);
 437       // We need to break here, or we could possibly update zemEn again
 438       break;
 439     }
 440     if (fillVz != kInvalidVtxZ) fHVertexZ->Fill(fillVz);
 441   }
 442   fZemEnergy->Fill(zemEn);
 443   fZdcEnergy->Fill(zdcEn);
 444   fCorrelationZemZdc->Fill(zemEn, zdcEn);
 445
 446   // --- Check if this is an outlier event ---------------------------
 447   if (CheckOutlier(ivtx, esd)) {
 448     fVertexZ = kInvalidVtxZ;
 449     return false;
 450   }
 451
 452   // --- Calculate the centrality ------------------------------------
 453   Float_t c1, c2, c3, c4, c5;
 454   Int_t runnumber = esd->GetRunNumber();
 455 #if 1 // New centrality determination by HHD
 456   if (runnumber < 137165 ) {
 457     // ref 137161
 458     c1 = 16992;
 459     c2 = 345;
 460     c3 = 2.23902e+02;
 461     c4 = 1.56667;
 462     c5 = 9.49434e-05;
 463   }
 464   else if (runnumber <= 137848 && runnumber >= 137230) {
 465     // ref 137722
 466     c1  = 15000;
 467     c2  = 295;
 468     c3  = 2.23660e+02;
 469     c4  = 1.56664;
 470     c5  = 8.99571e-05;
 471   }
 472   else if (runnumber <= 138275 && runnumber >= 138125) {
 473     // ref 137161, yes that's the correct run!
 474     c1 = 16992;
 475     c2 = 345;
 476     c3 = 2.23902e+02;
 477     c4 = 1.56667;
 478     c5 = 9.49434e-05;
 479   }
 480   else { // if (runnumber <= 139517 && runnumber >= 138358) {
 481     // Ref run 139172
 482     c1 = 16992;
 483     c2 = 345;
 484     c3 = 2.04789e+02;
 485     c4 = 1.56629;
 486     c5 = 1.02768e-04;
 487   }
 488 #else // Old centrality determination
 489   if (runnumber < 137722 && runnumber >= 137161 ) {
 490     c1 = 16992;
 491     c2 = 345;
 492     c3 = 2.23902e+02;
 493     c4 = 1.56667;
 494     c5 = 9.49434e-05;
 495   }
 496   else if (runnumber < 139172 && runnumber >= 137722) {
 497     c1  = 15000;
 498     c2  = 295;
 499     c3  = 2.23660e+02;
 500     c4 = 1.56664;
 501     c5 = 8.99571e-05;
 502   }
 503   else if (runnumber >= 139172) {
 504     c1 = 16992;
 505     c2 = 345;
 506     c3 = 2.04789e+02;
 507     c4 = 1.56629;
 508     c5 = 1.02768e-04;
 509   }
 510   else {
 511     c1 = 16992;
 512     c2 = 345;
 513     c3 = 2.04789e+02;
 514     c4 = 1.56629;
 515     c5 = 1.02768e-04;
 516   }
 517 #endif
 518   if (zemEn > c2) {
 519     Float_t slope   = (zdcEn + c1) / (zemEn - c2) + c3;
 520     Float_t zdcCent = (TMath::ATan(slope) - c4) / c5;
 521     if (zdcCent >= 0) fCent = zdcCent;
 522   }
 523
 524   fHCent->Fill(fCent);
 525
 526   return true;
 527 }
 528 //____________________________________________________________________
 529 //
 530 // EOF
 531 //