PWGLF/FORWARD/analysis2/AliDisplacedVertexSelection.cxx

   1 #include "AliDisplacedVertexSelection.h"
   2 #include <iostream>
   3 #include <TROOT.h>
   4 #include "AliESDEvent.h"
   5 #include "AliESDZDC.h"
   6 ClassImp(AliDisplacedVertexSelection)
   7 #if 0
   8 ; // This is for Emacs
   9 #endif
  10
  11 //____________________________________________________________________
  12 AliDisplacedVertexSelection::AliDisplacedVertexSelection()
  13   : TObject(),
  14     fVertexZ(9999),
  15     fCent(100)
  16 {
  17 }
  18 //____________________________________________________________________
  19 AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
  20   : TObject(o),
  21     fVertexZ(9999),
  22     fCent(100)
  23 {
  24 }
  25 //____________________________________________________________________
  26 AliDisplacedVertexSelection&
  27 AliDisplacedVertexSelection::operator=(const AliDisplacedVertexSelection& o)
  28 {
  29   if (&o == this) return *this;
  30   return *this;
  31 }
  32
  33 //____________________________________________________________________
  34 void
  35 AliDisplacedVertexSelection::Output(TList* /*l*/, const char* /* name*/) const
  36 {
  37 }
  38
  39 //____________________________________________________________________
  40 void
  41 AliDisplacedVertexSelection::Print(Option_t*) const
  42 {
  43 #if 0
  44   char ind[gROOT->GetDirLevel()+1];
  45   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  46   ind[gROOT->GetDirLevel()] = '\0';
  47   std::cout << std::boolalpha
  48             << std::noboolalpha << std::endl;
  49 #endif
  50 }
  51
  52 //____________________________________________________________________
  53 Bool_t
  54 AliDisplacedVertexSelection::Process(const AliESDEvent* esd)
  55 {
  56   fVertexZ = 9999; // Default vertex value
  57   fCent    = 100;  // Default centrality value
  58
  59   // Some constants
  60   const Float_t kZDCrefSum        = -66.5;
  61   const Float_t kZDCrefDelta      =  -2.10;
  62   const Float_t kZDCsigmaSum      =   3.25;
  63   const Float_t kZDCsigmaDelta    =   2.25;
  64   const Float_t kZDCsigmaSumSat   =   2.50;
  65   const Float_t kZDCsigmaDeltaSat =   1.35;
  66
  67   // Corrections for magnetic fields
  68   const Float_t kZEMcorrPlusPlus[21]   = {0.6840,0.7879,0.8722,
  69                                           0.9370,0.9837,1.0137,
  70                                           1.0292,1.0327,1.0271,
  71                                           1.0152,1.0000,0.9844,
  72                                           0.9714,0.9634,0.9626,
  73                                           0.9708,0.9891,1.0181,
  74                                           1.0574,1.1060,1.1617};
  75   const Float_t kZEMcorrMoinsMoins[21] = {0.7082,0.8012,0.8809,
  76                                           0.9447,0.9916,1.0220,
  77                                           1.0372,1.0395,1.0318,
  78                                           1.0174,1.0000,0.9830,
  79                                           0.9699,0.9635,0.9662,
  80                                           0.9794,1.0034,1.0371,
  81                                           1.0782,1.1224,1.1634};
  82
  83   // --- Get the ESD object ------------------------------------------
  84   AliESDZDC* esdZDC = esd->GetESDZDC();
  85   if (!esdZDC) {
  86     AliWarning("No ZDC ESD object!");
  87     return false;
  88   }
  89   Double_t currentL3   = esd->GetCurrentL3();
  90   Double_t currentDipo = esd->GetCurrentDip();
  91
  92   // --- ZDC and ZEM energy signal -----------------------------------
  93   Double_t zdcEn      = (esdZDC->GetZDCN1Energy()+
  94                          esdZDC->GetZDCP1Energy()+
  95                          esdZDC->GetZDCN2Energy()+
  96                          esdZDC->GetZDCP2Energy());
  97   Double_t zemEn      = (esdZDC->GetZDCEMEnergy(0)+
  98                          esdZDC->GetZDCEMEnergy(1))/8.;
  99
 100   // --- Calculate DeltaT and sumT -----------------------------------
 101   Double_t deltaTdc = 0;
 102   Double_t sumTdc   = 0;
 103
 104   for(Int_t i = 0; i < 4; ++i) {
 105     if(esdZDC->GetZDCTDCData(10,i) != 0) {
 106       Double_t  tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
 107                                     esdZDC->GetZDCTDCData(14,i));
 108       Double_t  tdcC       = esdZDC->GetZDCTDCCorrected(10,i);
 109       for(Int_t j = 0; j < 4; ++j) {
 110         if(esdZDC->GetZDCTDCData(12,j) != 0) {
 111           Double_t   tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
 112                                          esdZDC->GetZDCTDCData(14,j));
 113           Double_t   tdcA       = esdZDC->GetZDCTDCCorrected(12,j);
 114           if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
 115             deltaTdc = tdcC-tdcA;
 116             sumTdc   = tdcC+tdcA;
 117           }
 118           else {
 119             deltaTdc = tdcCnoCorr-tdcAnoCorr;
 120             sumTdc   = tdcCnoCorr+tdcAnoCorr;
 121           }
 122         }
 123       }
 124     }
 125   }
 126
 127   // --- Find the vertex ---------------------------------------------
 128   if(deltaTdc!=0. || sumTdc!=0.) {
 129     for (Int_t k = -10; k <= 10; ++k) {
 130       Float_t zsat  = 2.55F * k;
 131       Float_t delta = (k == 0 ? kZDCsigmaDelta : kZDCsigmaDeltaSat);
 132       Float_t sum   = (k == 0 ? kZDCsigmaSum   : kZDCsigmaSumSat);
 133       Float_t dT    = deltaTdc - kZDCrefDelta - zsat;
 134       Float_t sT    = sumTdc  - kZDCrefSum - zsat;
 135       Float_t check = dT * dT / delta / delta + sT * sT / sum  / sum;
 136       if (check > 1.0 || k == 0) continue;
 137
 138       // Set the vertex
 139       fVertexZ = 37.5 * k;
 140
 141       // Correct zem energy
 142       if(currentDipo>0 && currentL3>0) zemEn /= kZEMcorrPlusPlus[k+10];
 143       if(currentDipo<0 && currentL3<0) zemEn /= kZEMcorrMoinsMoins[k+10];
 144     }
 145   }
 146
 147   // --- Calculate the centrality ------------------------------------
 148   Float_t c1, c2, c3, c4, c5;
 149   Int_t runnumber = esd->GetRunNumber();
 150   if (runnumber < 137722 && runnumber >= 137161 ) {
 151     c1 = 16992;
 152     c2 = 345;
 153     c3 = 2.23902e+02;
 154     c4 = 1.56667;
 155     c5 = 9.49434e-05;
 156   }
 157   else if (runnumber < 139172 && runnumber >= 137722) {
 158     c1  = 15000;
 159     c2  = 295;
 160     c3  = 2.23660e+02;
 161     c4 = 1.56664;
 162     c5 = 8.99571e-05;
 163   }
 164   else if (runnumber >= 139172) {
 165     c1 = 16992;
 166     c2 = 345;
 167     c3 = 2.04789e+02;
 168     c4 = 1.56629;
 169     c5 = 1.02768e-04;
 170   }
 171   else {
 172     c1 = 16992;
 173     c2 = 345;
 174     c3 = 2.04789e+02;
 175     c4 = 1.56629;
 176     c5 = 1.02768e-04;
 177   }
 178
 179   if (zemEn > c2) {
 180     Float_t slope   = (zdcEn + c1) / (zemEn - c2) + c3;
 181     Float_t zdcCent = (TMath::ATan(slope) - c4) / c5;
 182     if (zdcCent >= 0) fCent = zdcCent;
 183   }
 184
 185   return true;
 186 }
 187
 188 //____________________________________________________________________
 189 Double_t
 190 AliDisplacedVertexSelection::CheckDisplacedVertex(const AliESDEvent* esd) const {
 191   // Code taken directly from M.Guilbaud.
 192   AliWarning("This method is deprecated");
 193   Double_t zvtx = 9999.;
 194
 195   Float_t kZDCrefSum        =-66.5;
 196   Float_t kZDCrefDelta      =-2.10;
 197   Float_t kZDCsigmaSum      = 3.25;
 198   Float_t kZDCsigmaDelta    = 2.25;
 199   Float_t kZDCsigmaSumSat   = 2.50;
 200   Float_t kZDCsigmaDeltaSat = 1.35;
 201
 202   AliESDZDC* esdZDC = esd -> GetESDZDC();
 203
 204   /* Double_t zdcNCEnergy  = esdZDC->GetZDCN1Energy();
 205      Double_t fZpcEnergy  = esdZDC->GetZDCP1Energy();
 206      Double_t fZnaEnergy  = esdZDC->GetZDCN2Energy();
 207      Double_t fZpaEnergy  = esdZDC->GetZDCP2Energy();
 208      Double_t fZem1Energy = esdZDC->GetZDCEMEnergy(0);
 209      Double_t fZem2Energy = esdZDC->GetZDCEMEnergy(1);*/
 210
 211   //Double_t fzdcEn      = (esdZDC->GetZDCN1Energy()+esdZDC->GetZDCP1Energy()+esdZDC->GetZDCN2Energy()+esdZDC->GetZDCP2Energy());
 212   //Double_t fzemEn      = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.;
 213
 214
 215
 216   ///////////////////
 217   //Event Selection//
 218   ///////////////////
 219   Double_t deltaTdc = 0;
 220   Double_t sumTdc   = 0;
 221
 222   for(Int_t i = 0; i < 4; ++i) {
 223     if(esdZDC->GetZDCTDCData(10,i) != 0) {
 224       Double_t  tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
 225                                     esdZDC->GetZDCTDCData(14,i));
 226       Double_t  tdcC       = esdZDC->GetZDCTDCCorrected(10,i);
 227       for(Int_t j = 0; j < 4; ++j) {
 228         if(esdZDC->GetZDCTDCData(12,j) != 0) {
 229           Double_t   tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
 230                                          esdZDC->GetZDCTDCData(14,j));
 231           Double_t   tdcA       = esdZDC->GetZDCTDCCorrected(12,j);
 232           if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
 233             deltaTdc = tdcC-tdcA;
 234             sumTdc   = tdcC+tdcA;
 235           }
 236           else {
 237             deltaTdc = tdcCnoCorr-tdcAnoCorr;
 238             sumTdc   = tdcCnoCorr+tdcAnoCorr;
 239           }
 240         }
 241       }
 242     }
 243   }
 244
 245   ///////////////////////
 246   //Global Event Filter//
 247   ///////////////////////
 248   Bool_t zdcAccSat[21];
 249
 250   // Bool_t zdcAccSatRunClass[21];
 251   for(Int_t k = -10; k <= 10; k++) zdcAccSat[k+10] = kFALSE;
 252
 253
 254   if(deltaTdc!=0. || sumTdc!=0.) {
 255     for (Int_t k = -10; k <= 10; ++k) {
 256       Float_t zsat = 2.55F * k;
 257
 258       if(k==0) {
 259         if(((deltaTdc-kZDCrefDelta)*(deltaTdc-kZDCrefDelta)/
 260             (kZDCsigmaDelta*kZDCsigmaDelta)+
 261             (sumTdc  -kZDCrefSum  ) * (sumTdc  -kZDCrefSum  ) /
 262             (kZDCsigmaSum  *kZDCsigmaSum))<=1.0) {
 263           zdcAccSat[k+10] = kTRUE;
 264         }
 265       }
 266       else {
 267         if(((deltaTdc-kZDCrefDelta-zsat)*(deltaTdc-kZDCrefDelta-zsat)/
 268             (kZDCsigmaDeltaSat*kZDCsigmaDeltaSat)+
 269             (sumTdc  -kZDCrefSum  -zsat)*(sumTdc  -kZDCrefSum  -zsat)/
 270             (kZDCsigmaSumSat  *kZDCsigmaSumSat  ))<=1.0) {
 271           zdcAccSat[k+10] = kTRUE;
 272         }
 273       }
 274     }
 275   }
 276
 277   for(Int_t k=-10;k<=10;++k) {
 278     if(zdcAccSat[k+10] && k!=0 ) {
 279       //std::cout<<"Displaced vertex at "<<37.5*(Float_t)k<<" cm"<<std::endl;
 280       zvtx = 37.5*(Float_t)k;
 281     }
 282   }
 283
 284   return zvtx;
 285
 286 }
 287 //____________________________________________________________________
 288 Double_t
 289 AliDisplacedVertexSelection::CalculateDisplacedVertexCent(const AliESDEvent* esd) const
 290 {
 291   Float_t kZDCrefSum        =-66.5;
 292   Float_t kZDCrefDelta      =-2.10;
 293   Float_t kZDCsigmaSum      = 3.25;
 294   Float_t kZDCsigmaDelta    = 2.25;
 295   Float_t kZDCsigmaSumSat   = 2.50;
 296   Float_t kZDCsigmaDeltaSat = 1.35;
 297
 298   AliESDZDC* esdZDC = esd -> GetESDZDC();
 299
 300   Int_t runnumber = esd->GetRunNumber();
 301   Double_t fcurrentL3 = esd->GetCurrentL3();
 302   Double_t fcurrentDipo = esd->GetCurrentDip();
 303   Double_t cent = -1;
 304
 305   Double_t fzdcEn      = (esdZDC->GetZDCN1Energy()+
 306                           esdZDC->GetZDCP1Energy()+
 307                           esdZDC->GetZDCN2Energy()+
 308                           esdZDC->GetZDCP2Energy());
 309   Double_t fzemEn      = (esdZDC->GetZDCEMEnergy(0)+
 310                           esdZDC->GetZDCEMEnergy(1))/8.;
 311
 312
 313   ///////////////////
 314   //Event Selection//
 315   ///////////////////
 316   Double_t deltaTdc = 0;
 317   Double_t sumTdc   = 0;
 318
 319   for(Int_t i = 0; i < 4; ++i) {
 320     if(esdZDC->GetZDCTDCData(10,i) != 0) {
 321       Double_t  tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
 322                                     esdZDC->GetZDCTDCData(14,i));
 323       Double_t  tdcC       = esdZDC->GetZDCTDCCorrected(10,i);
 324       for(Int_t j = 0; j < 4; ++j) {
 325         if(esdZDC->GetZDCTDCData(12,j) != 0) {
 326           Double_t   tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
 327                                          esdZDC->GetZDCTDCData(14,j));
 328           Double_t   tdcA       = esdZDC->GetZDCTDCCorrected(12,j);
 329           if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
 330             deltaTdc = tdcC-tdcA;
 331             sumTdc   = tdcC+tdcA;
 332           }
 333           else {
 334             deltaTdc = tdcCnoCorr-tdcAnoCorr;
 335             sumTdc   = tdcCnoCorr+tdcAnoCorr;
 336           }
 337         }
 338       }
 339     }
 340   }
 341
 342   ///////////////////////
 343   //Global Event Filter//
 344   ///////////////////////
 345   Bool_t zdcAccSat[21];
 346   // Bool_t zdcAccSatRunClass[21];
 347   for(Int_t k = -10; k <= 10; k++)  zdcAccSat[k+10]         = kFALSE;
 348
 349   if(deltaTdc!=0. || sumTdc!=0.) {
 350     for(Int_t k = -10; k <= 10; ++k) {
 351       Float_t zsat = 2.55*(Float_t)k;
 352
 353       if(k==0) {
 354         if(((deltaTdc-kZDCrefDelta)*(deltaTdc-kZDCrefDelta)/
 355             (kZDCsigmaDelta*kZDCsigmaDelta)+
 356             (sumTdc  -kZDCrefSum  )*(sumTdc  -kZDCrefSum  )/
 357             (kZDCsigmaSum  *kZDCsigmaSum))<=1.0) {
 358           zdcAccSat[k+10] = kTRUE;
 359         }
 360       }
 361       else {
 362         if(((deltaTdc-kZDCrefDelta-zsat)*(deltaTdc-kZDCrefDelta-zsat)/
 363             (kZDCsigmaDeltaSat*kZDCsigmaDeltaSat)+
 364             (sumTdc  -kZDCrefSum  -zsat)*(sumTdc  -kZDCrefSum  -zsat)/
 365             (kZDCsigmaSumSat  *kZDCsigmaSumSat  ))<=1.0) {
 366           zdcAccSat[k+10] = kTRUE;
 367         }
 368
 369       }
 370
 371     }
 372   }
 373
 374   Float_t kZEMcorrPlusPlus[21]   = {0.6840,0.7879,0.8722,0.9370,0.9837,1.0137,
 375                                     1.0292,1.0327,1.0271,1.0152,1.0000,0.9844,
 376                                     0.9714,0.9634,0.9626,0.9708,0.9891,1.0181,
 377                                     1.0574,1.1060,1.1617};
 378   Float_t kZEMcorrMoinsMoins[21] = {0.7082,0.8012,0.8809,0.9447,0.9916,1.0220,
 379                                     1.0372,1.0395,1.0318,1.0174,1.0000,0.9830,
 380                                     0.9699,0.9635,0.9662,0.9794,1.0034,1.0371,
 381                                     1.0782,1.1224,1.1634};
 382   ///////////////////
 383   //ZemEn correction//
 384   ////////////////////
 385
 386   for(Int_t k = -10; k <= 10; ++k) {
 387     if(zdcAccSat[k+10]) {
 388       //std::cout<<"Vertex at "<<37.5*(Float_t)k<<"cm from cent"<<std::endl;
 389       if(fcurrentDipo>0 && fcurrentL3>0) {
 390         fzemEn /= kZEMcorrPlusPlus[k+10];
 391       }
 392       if(fcurrentDipo<0 && fcurrentL3<0) {
 393         fzemEn /= kZEMcorrMoinsMoins[k+10];
 394       }
 395     }
 396   }
 397
 398   ////////////////////////
 399   //Centrality selection//
 400   ////////////////////////
 401   Double_t fzdcPercentile = -1;
 402   Float_t slope;
 403
 404
 405   if(runnumber < 137722 && runnumber >= 137161 ) {
 406     if(fzemEn > 345.) {
 407       slope = (fzdcEn + 16992.)/(fzemEn - 345.);
 408       slope += 2.23902e+02;
 409       fzdcPercentile = (TMath::ATan(slope) - 1.56667)/9.49434e-05;
 410       if (fzdcPercentile<0) fzdcPercentile = 0;
 411     }
 412     else fzdcPercentile = 100;                }
 413   else if(runnumber < 139172 && runnumber >= 137722) {
 414     if(fzemEn > 295.) {
 415       slope = (fzdcEn + 15000.)/(fzemEn - 295.);
 416       slope += 2.23660e+02;
 417       fzdcPercentile = (TMath::ATan(slope) - 1.56664)/8.99571e-05;
 418       if (fzdcPercentile<0) fzdcPercentile = 0;
 419     }
 420     else fzdcPercentile = 100;
 421   }
 422   else if(runnumber >= 139172) {
 423     if(fzemEn > 345.) {
 424       slope = (fzdcEn + 16992.)/(fzemEn - 345.);
 425       slope += 2.04789e+02;
 426       fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
 427       if (fzdcPercentile<0) fzdcPercentile = 0;
 428     }
 429     else fzdcPercentile = 100;
 430   }
 431   else {
 432     if(fzemEn > 345.) {
 433       slope = (fzdcEn + 16992.)/(fzemEn - 345.);
 434       slope += 2.04789e+02;
 435       fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
 436       if(fzdcPercentile<0) fzdcPercentile = 0;
 437     }
 438     else fzdcPercentile = 100;
 439   }
 440
 441   if(fzdcPercentile > 0 && fzdcPercentile <100) {
 442     //std::cout<<fzdcPercentile<<std::endl;
 443     cent = fzdcPercentile;
 444   }
 445   return cent;
 446 }
 447 //____________________________________________________________________
 448 //
 449 // EOF
 450 //