PWGLF/FORWARD/analysis2/AliDisplacedVertexSelection.cxx

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