AliPhysics  f05a842 (f05a842)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliDisplacedVertexSelection.cxx
Go to the documentation of this file.
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"
15 #if 0
16 ; // This is for Emacs
17 #endif
18 
19 //____________________________________________________________________
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 //____________________________________________________________________
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 //____________________________________________________________________
54 {
55  if (&o == this) return *this;
56 
57  fDeltaTdc = o.fDeltaTdc;
58  fSumTdc = o.fSumTdc;
63  fMC = o.fMC;
64 
65  return *this;
66 }
67 
68 //____________________________________________________________________
69 void
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
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
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
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
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)) {
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 //
Bool_t Process(const AliESDEvent *esd)
double Double_t
Definition: External.C:58
int Int_t
Definition: External.C:63
float Float_t
Definition: External.C:68
Bool_t ProcessMC(const AliMCEvent *mcevent)
Definition: External.C:228
Definition: External.C:212
Float_t GetZemCorr(Int_t k, Bool_t minusminus) const
Selection of events from satellite interactions.
ClassImp(AliDisplacedVertexSelection) AliDisplacedVertexSelection
AliDisplacedVertexSelection & operator=(const AliDisplacedVertexSelection &o)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void SetupForData(TList *l, const char *name=0, Bool_t mc=false)
Bool_t CheckOutlier(Int_t ivtx, const AliESDEvent *esd) const
void Print(Option_t *option="") const