3 #include "AliMCEvent.h"
4 #include "AliMultiplicity.h"
7 #include "AliGenEventHeader.h"
12 #include "AliESDEvent.h"
13 #include "AliESDZDC.h"
26 fCorrelationZemZdc(0),
27 fCorrelationSumDelta(0),
28 fVertexZ(kInvalidVtxZ),
38 fDeltaTdc(o.fDeltaTdc),
40 fZdcEnergy(o.fZdcEnergy),
41 fZemEnergy(o.fZemEnergy),
42 fCorrelationZemZdc(o.fCorrelationZemZdc),
43 fCorrelationSumDelta(o.fCorrelationSumDelta),
44 fVertexZ(kInvalidVtxZ),
55 if (&o ==
this)
return *
this;
77 out->SetName(
"displacedVertex");
86 2*
kMaxK+1, vzMin, vzMax);
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");
100 fHCent->SetFillColor(kBlue+1);
101 fHCent->SetFillStyle(3001);
104 Int_t nDeltaT = 2000;
108 nDeltaT,minDeltaT,maxDeltaT);
113 Int_t nSumT = nDeltaT;
116 fSumTdc =
new TH1D(
"SumTdc",
"#sumTDC",nSumT,minSumT,maxSumT);
124 fZdcEnergy =
new TH1D(
"ZDCEnergy",
"E_{ZDC}",nZdcE,minZdcE,maxZdcE);
132 fZemEnergy =
new TH1D(
"ZEMEnergy",
"E_{ZEM}",nZemE,minZemE,maxZemE);
138 nZemE, minZemE, maxZemE,
139 nZdcE, minZdcE, maxZdcE);
146 nSumT, minSumT, maxSumT,
147 nDeltaT, minDeltaT, maxDeltaT);
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;
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};
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 };
209 return minusminus ? kMoinsMoinsMC[i] : kPlusPlusMC[i];
210 return minusminus ? kMoinsMoins[i] : kPlusPlus[i];
219 if (
fMC)
return false;
220 if (ivtx <= 4)
return false;
223 AliESDVZERO* esdV0 = esd->GetVZEROData();
224 Double_t nClusters = esd->GetMultiplicity()->GetNumberOfITSClusters(1);
228 const Double_t slp[8][21] = {{0.000,0.000,0.000,0.000,0.000,0.000,0.000,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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,
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}};
279 for (
int iring = 4; iring < 8; iring++) {
281 for (
int iCh=0; iCh<8; iCh++) {
282 Int_t idx = iCh+iring*8;
283 multRing += esdV0->GetMultiplicity(idx);
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);
291 Double_t slpDOWN = slp[iring][ivtx]*(1.-downFac);
292 Double_t constant = cst[iring][ivtx];
294 if ((slpDOWN * nClusters + constant) > multRing ||
295 (slpUP * nClusters + constant) < multRing ) {
307 if (!
fMC)
return true;
308 if (!mcevent)
return false;
310 AliMCEvent*
event =
const_cast<AliMCEvent*
>(mcevent);
311 AliHeader* header =
event->Header();
312 AliGenEventHeader* genHeader = header->GenEventHeader();
314 genHeader->PrimaryVertex(vertex);
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.};
324 for(
Int_t iVtx=0; iVtx<16; ++iVtx) {
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))) {
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;
358 AliESDZDC* esdZDC = esd->GetESDZDC();
360 AliWarning(
"No ZDC ESD object!");
363 Double_t currentL3 = esd->GetCurrentL3();
364 Double_t currentDipo = esd->GetCurrentDip();
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.;
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));
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;
403 deltaTdc = tdcCnoCorr-tdcAnoCorr;
404 sumTdc = tdcCnoCorr+tdcAnoCorr;
416 if(deltaTdc!=0. || sumTdc!=0.) {
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;
435 if(currentDipo>0 && currentL3>0) zemEn /=
GetZemCorr(k,
false);
436 if(currentDipo<0 && currentL3<0) zemEn /=
GetZemCorr(k,
true);
454 Int_t runnumber = esd->GetRunNumber();
455 #if 1 // New centrality determination by HHD
456 if (runnumber < 137165 ) {
464 else if (runnumber <= 137848 && runnumber >= 137230) {
472 else if (runnumber <= 138275 && runnumber >= 138125) {
488 #else // Old centrality determination
489 if (runnumber < 137722 && runnumber >= 137161 ) {
496 else if (runnumber < 139172 && runnumber >= 137722) {
503 else if (runnumber >= 139172) {
519 Float_t slope = (zdcEn + c1) / (zemEn - c2) + c3;
520 Float_t zdcCent = (TMath::ATan(slope) - c4) / c5;
521 if (zdcCent >= 0)
fCent = zdcCent;
Bool_t Process(const AliESDEvent *esd)
TH2D * fCorrelationZemZdc
Bool_t ProcessMC(const AliMCEvent *mcevent)
AliDisplacedVertexSelection()
Float_t GetZemCorr(Int_t k, Bool_t minusminus) const
Selection of events from satellite interactions.
ClassImp(AliDisplacedVertexSelection) AliDisplacedVertexSelection
AliDisplacedVertexSelection & operator=(const AliDisplacedVertexSelection &o)
void SetupForData(TList *l, const char *name=0, Bool_t mc=false)
TH2D * fCorrelationSumDelta
Bool_t CheckOutlier(Int_t ivtx, const AliESDEvent *esd) const
void Print(Option_t *option="") const