AliPhysics  2c6b7ad (2c6b7ad)
AliEmcalCorrectionCellEmulateCrosstalk.cxx
Go to the documentation of this file.
1 // AliEmcalCorrectionCellEmulateCrosstalk
2 //
3 
4 #include <map>
5 #include <vector>
6 #include <string>
7 
8 #include <TObjArray.h>
9 #include <TFile.h>
10 
11 #include <AliEMCALGeometry.h>
12 #include <AliEMCALRecoUtils.h>
13 #include <AliAODEvent.h>
14 
16 
20 
21 // Actually registers the class with the base class
23 
28  AliEmcalCorrectionComponent("AliEmcalCorrectionCellEmulateCrosstalk"),
29  fCellEnergyDistBefore(0),
30  fCellEnergyDistAfter(0),
31  fTCardCorrClusEnerConserv(kFALSE),
32  fRandom(0),
33  fRandomizeTCard(kTRUE),
34  fTCardCorrMinAmp(0.01),
35  fTCardCorrMinInduced(0),
36  fTCardCorrMaxInducedELeak(0),
37  fTCardCorrMaxInduced(100),
38  fPrintOnce(kFALSE),
39  fAODCellsTmp(0x0)
40 {
41  for(Int_t i = 0; i < fgkNsm; i++)
42  {
46 
47  for(Int_t j = 0; j < 4 ; j++)
48  {
49  fTCardCorrInduceEner [j][i] = 0 ;
50  fTCardCorrInduceEnerFrac [j][i] = 0 ;
51  fTCardCorrInduceEnerFracP1 [j][i] = 0 ;
53  }
54  }
55 
56  ResetArrays();
57 }
58 
63 {
64  fAODCellsTmp->DeleteContainer();
65  delete fAODCellsTmp;
66 }
67 
72 {
73  // Initialization
75 
76  AliWarning("Init EMCAL crosstalk emulation");
77 
78  GetProperty("conservEnergy", fTCardCorrClusEnerConserv);
79  GetProperty("randomizeTCardInducedEnergy", fRandomizeTCard);
80  GetProperty("inducedTCardMinimumCellEnergy", fTCardCorrMinAmp);
81  GetProperty("inducedTCardMaximum", fTCardCorrMaxInduced);
82  GetProperty("inducedTCardMaximumELeak", fTCardCorrMaxInducedELeak);
83  GetProperty("inducedTCardMinimum", fTCardCorrMinInduced);
84 
85  // Handle array initialization
86  // For the format of these values, see the "default" yaml configuration file
87  const std::map <std::string, Float_t (*)[fgkNsm]> properties2D = {
88  {"inducedEnergyLossConstant", fTCardCorrInduceEner},
89  {"inducedEnergyLossFraction", fTCardCorrInduceEnerFrac},
90  {"inducedEnergyLossFractionP1", fTCardCorrInduceEnerFracP1},
91  {"inducedEnergyLossFractionWidth", fTCardCorrInduceEnerFracWidth}
92  };
93  const std::map <std::string, Float_t *> properties1D = {
94  {"inducedEnergyLossMinimumFraction", fTCardCorrInduceEnerFracMin},
95  {"inducedEnergyLossMaximumFraction", fTCardCorrInduceEnerFracMax},
96  {"inducedEnergyLossProbability", fTCardCorrInduceEnerProb}
97  };
98  RetrieveAndSetProperties(properties2D);
99  RetrieveAndSetProperties(properties1D);
100 
101  if (!fRecoUtils) {
103  }
104 
105  GetProperty("printConfiguration", fPrintOnce);
106  if (fPrintOnce) {
107  PrintTCardParam();
108  }
109 
110  return kTRUE;
111 }
112 
121 void AliEmcalCorrectionCellEmulateCrosstalk::SetProperty(Float_t val[][fgkNsm], std::vector<double> & property, unsigned int iSM, const std::string & name)
122 {
123  for (unsigned int iProperty = 0; iProperty < property.size(); iProperty++) {
124  val[iProperty][iSM] = property.at(iProperty);
125  }
126 }
127 
136 void AliEmcalCorrectionCellEmulateCrosstalk::SetProperty(Float_t val[fgkNsm], std::vector<double> & property, unsigned int iSM, const std::string & name)
137 {
138  if (property.size() != 1) {
139  AliErrorStream() << "Trying to set single value for property " << name << ", but given " << property.size() << " values. Please check your configuration!\n";
140  }
141  else {
142  val[iSM] = property.at(0);
143  }
144 }
145 
150 {
152 
153  if (fCreateHisto){
154  fCellEnergyDistBefore = new TH1F("hCellEnergyDistBefore","hCellEnergyDistBefore;E_{cell} (GeV)",1000,0,10);
156  fCellEnergyDistAfter = new TH1F("hCellEnergyDistAfter","hCellEnergyDistAfter;E_{cell} (GeV)",1000,0,10);
158  }
159 }
160 
165 {
167 
168  if (!fEventManager.InputEvent()) {
169  AliError("Event ptr = 0, returning");
170  return kFALSE;
171  }
172 
173  // START PROCESSING ---------------------------------------------------------
174  // Test if cells present
175  if (fCaloCells->GetNumberOfCells()<=0)
176  {
177  AliDebug(2, Form("Number of EMCAL cells = %d, returning", fCaloCells->GetNumberOfCells()));
178  return kFALSE;
179  }
180 
181  if(fCreateHisto)
182  FillCellQA(fCellEnergyDistBefore); // "before" QA
183 
184  // CELL CROSSTALK EMULATION -------------------------------------------------------
185 
186  // Compute the induced cell energies by T-Card correlation emulation, ONLY MC
188 
189  // Add to existing cells the found induced energies in MakeCellTCardCorrelation() if new signal is larger than 10 MeV.
191 
192  // Add new cells with found induced energies in MakeCellTCardCorrelation() if new signal is larger than 10 MeV.
194 
195  // -------------------------------------------------------
196 
197  if(fCreateHisto)
198  FillCellQA(fCellEnergyDistAfter); // "after" QA
199 
200  ResetArrays();
201 
202  return kTRUE;
203 }
204 
210 {
211  Int_t id = -1;
212  Float_t amp = -1;
213 
214  if (fPrintOnce) {
215  PrintTCardParam();
216  fPrintOnce = 0;
217  }
218 
219  Int_t nCells = fCaloCells->GetNumberOfCells();
220 
221  // Loop on all cells with signal
222  for (Int_t icell = 0; icell < nCells; icell++)
223  {
224  id = fCaloCells->GetCellNumber(icell);
225  amp = fCaloCells->GetAmplitude (icell); // fCaloCells->GetCellAmplitude(id);
226 
227  if ( amp <= fTCardCorrMinAmp ) continue ;
228 
229  //
230  // First get the SM, col-row of this tower
231  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
232  fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
233  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
234 
235  //
236  // Determine randomly if we want to create a correlation for this cell,
237  // depending the SM number of the cell
238  if ( fTCardCorrInduceEnerProb[imod] < 1 )
239  {
240  Float_t rand = fRandom.Uniform(0, 1);
241 
242  if ( rand > fTCardCorrInduceEnerProb[imod] ) continue;
243  }
244 
245  AliDebug(1,Form("Reference cell absId %d, iEta %d, iPhi %d, amp %2.3f",id,ieta,iphi,amp));
246 
247  //
248  // Get the absId of the cells in the cross and same T-Card
249  Int_t absIDup = -1;
250  Int_t absIDdo = -1;
251  Int_t absIDlr = -1;
252  Int_t absIDuplr = -1;
253  Int_t absIDdolr = -1;
254 
255  Int_t absIDup2 = -1;
256  Int_t absIDup2lr = -1;
257  Int_t absIDdo2 = -1;
258  Int_t absIDdo2lr = -1;
259 
260  // Only 2 columns in the T-Card, +1 for even and -1 for odd with respect reference cell
261  Int_t colShift = 0;
262  if ( (ieta%2) && ieta <= AliEMCALGeoParams::fgkEMCALCols-1 ) colShift = -1;
263  if ( !(ieta%2) && ieta >= 0 ) colShift = +1;
264 
265  absIDlr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+colShift);
266 
267  // Check up / down cells from reference cell not out of SM and in same T-Card
268  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1 )
269  {
270  absIDup = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
271  absIDuplr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta+colShift);
272  }
273 
274  if ( iphi > 0 )
275  {
276  absIDdo = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
277  absIDdolr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta+colShift);
278  }
279 
280  // Check 2 up / 2 down cells from reference cell not out of SM
281  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-2 )
282  {
283  absIDup2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta);
284  absIDup2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta+colShift);
285  }
286 
287  if ( iphi > 1 )
288  {
289  absIDdo2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta);
290  absIDdo2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta+colShift);
291  }
292 
293  // In same T-Card?
294  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+1)/8) ) { absIDup = -1 ; absIDuplr = -1 ; }
295  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-1)/8) ) { absIDdo = -1 ; absIDdolr = -1 ; }
296  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+2)/8) ) { absIDup2 = -1 ; absIDup2lr = -1 ; }
297  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-2)/8) ) { absIDdo2 = -1 ; absIDdo2lr = -1 ; }
298 
299  // Calculate induced energy to T-Card cells
300  //
301  AliDebug(1,Form("cell up %d:" ,absIDup));
302  CalculateInducedEnergyInTCardCell(absIDup , id, imod, amp, 0);
303  AliDebug(1,Form("cell down %d:",absIDdo));
304  CalculateInducedEnergyInTCardCell(absIDdo , id, imod, amp, 0);
305 
306  AliDebug(1,Form("cell up left-right %d:" ,absIDuplr));
307  CalculateInducedEnergyInTCardCell(absIDuplr , id, imod, amp, 1);
308  AliDebug(1,Form("cell down left-right %d:",absIDdolr));
309  CalculateInducedEnergyInTCardCell(absIDdolr , id, imod, amp, 1);
310 
311  AliDebug(1,Form("cell left-right %d:",absIDlr));
312  CalculateInducedEnergyInTCardCell(absIDlr , id, imod, amp, 2);
313 
314  AliDebug(1,Form("cell up 2nd row %d:" ,absIDup2));
315  CalculateInducedEnergyInTCardCell(absIDup2 , id, imod, amp, 3);
316  AliDebug(1,Form("cell down 2nd row %d:",absIDdo2));
317  CalculateInducedEnergyInTCardCell(absIDdo2 , id, imod, amp, 3);
318 
319  AliDebug(1,Form("cell up left-right 2nd row %d:" ,absIDup2lr));
320  CalculateInducedEnergyInTCardCell(absIDup2lr, id, imod, amp, 3);
321  AliDebug(1,Form("cell down left-right 2nd row %d:",absIDdo2lr));
322  CalculateInducedEnergyInTCardCell(absIDdo2lr, id, imod, amp, 3);
323 
324  } // cell loop
325 
326 }
327 
328 //_______________________________________________________
337 //_______________________________________________________
339 (Int_t absId, Int_t absIdRef, Int_t sm, Float_t ampRef, Int_t cellCase)
340 {
341  // Check that the cell exists
342  if ( !AcceptCell(absId) ) return ;
343 
344  // Get the fraction
345  Float_t frac = fTCardCorrInduceEnerFrac[cellCase][sm] + ampRef * fTCardCorrInduceEnerFracP1[cellCase][sm];
346 
347  // Use an absolute minimum and maximum fraction if calculated one is out of range
348  if ( frac < fTCardCorrInduceEnerFracMin[sm] ) frac = fTCardCorrInduceEnerFracMin[sm];
349  if ( frac > fTCardCorrInduceEnerFracMax[sm] ) frac = fTCardCorrInduceEnerFracMax[sm];
350 
351  AliDebug(1,Form("\t fraction %2.3f",frac));
352 
353  // Randomize the induced fraction, if requested
354  if ( fRandomizeTCard )
355  {
356  frac = fRandom.Gaus(frac, fTCardCorrInduceEnerFracWidth[cellCase][sm]);
357 
358  AliDebug(1,Form("\t randomized fraction %2.3f",frac));
359  }
360 
361  // If too small or negative, do nothing else
362  if ( frac < 0.0001 ) return;
363 
364  // Calculate induced energy
365  Float_t inducedE = fTCardCorrInduceEner[cellCase][sm] + ampRef * frac;
366 
367  // Check if we induce too much energy, in such case use a constant value
368  if ( fTCardCorrMaxInduced < inducedE ) inducedE = fTCardCorrMaxInduced;
369 
370  AliDebug(1,Form("\t induced E %2.3f",inducedE));
371 
372  // Add the induced energy, check if cell existed
373  // Check that the induced+amp is large enough to avoid extra linearity effects
374  // typically of the order of the clusterization cell energy cut
375  // But if it is below 1 ADC, typically 10 MeV, also do it, to match Beam test linearity
376  Float_t amp = fCaloCells->GetCellAmplitude(absId) ;
377  if ( (amp+inducedE) > fTCardCorrMinInduced || inducedE < fTCardCorrMaxInducedELeak )
378  {
379  fTCardCorrCellsEner[absId] += inducedE;
380 
381  // If original energy of cell was null, create new one
382  if ( amp < 0.01 ) fTCardCorrCellsNew[absId] = kTRUE;
383  }
384  else return ;
385 
386  // Subtract the added energy to main cell, if energy conservation is requested
388  fTCardCorrCellsEner[absIdRef] -= inducedE;
389 }
390 
391 
398 {
399  Short_t absId = -1;
400  Double_t amp = -1;
401  Double_t time = -1;
402  Int_t mclabel = -1;
403  Double_t efrac = 0.;
404 
405  // Add the induced energy to the cells and copy them into a new temporal container
406  // used in AddInducedEnergiesToNewCells() to refill the default cells list fCaloCells
407  // Create the data member only once. Done here, not sure where to do this properly in the framework.
408  //
409  if ( !fAODCellsTmp )
410  fAODCellsTmp = new AliAODCaloCells("tmpCaloCellsAOD","tmpCaloCellsAOD",AliAODCaloCells::kEMCALCell);
411 
412  Int_t nCells = fCaloCells->GetNumberOfCells();
413  fAODCellsTmp->CreateContainer(nCells);
414 
415  for (Int_t icell = 0; icell < nCells; icell++)
416  {
417  // Get cell
418  fCaloCells->GetCell(icell, absId, amp, time, mclabel, efrac);
419 
420  amp+=fTCardCorrCellsEner[absId];
421 
422  // Set new amplitude in new temporal container
423  fAODCellsTmp->SetCell(icell, absId, amp, time, mclabel, efrac);
424  }
425 }
426 
431 {
432  // Count how many new cells
433  //
434  Int_t nCells = fCaloCells->GetNumberOfCells();
435  Int_t nCellsNew = 0;
436  for(Int_t j = 0; j < fgkNEMCalCells; j++)
437  {
438  // Newly created?
439  if ( !fTCardCorrCellsNew[j] ) continue;
440 
441  // Accept only if at least 10 MeV
442  if ( fTCardCorrCellsEner[j] < 0.01 ) continue;
443 
444  nCellsNew++;
445  }
446 
447  // Delete default cells container and create new one with larger arrays
448  // to contain new cells.
449  //
450  fCaloCells->DeleteContainer();
451  fCaloCells->CreateContainer(nCells+nCellsNew);
452 
453  // Recover the original input from fAODCellsTmp filled in AddInducedEnergiesToExistingCells()
454  //
455  Short_t absId = -1;
456  Float_t amp = -1;
457  Double_t time = 0;
458  Int_t mclabel = -1;
459  Double_t efrac = -1;
460  Bool_t highgain = 1;
461  Int_t cellNumber = -1;
462 
463  for(Int_t j = 0; j < nCells; j++)
464  {
465  absId = fAODCellsTmp->GetCellNumber(j);
466  amp = fAODCellsTmp->GetAmplitude(j);
467  time = fAODCellsTmp->GetTime(j);
468  mclabel = fAODCellsTmp->GetMCLabel(j);
469  efrac = fAODCellsTmp->GetEFraction(j);
470  highgain = fAODCellsTmp->GetHighGain(j);
471  fCaloCells->SetCell(j, absId, amp, time, mclabel, efrac, highgain);
472  }
473 
474  // Add the new cells
475  //
476  for(Int_t j = 0; j < fgkNEMCalCells; j++)
477  {
478  // Newly created?
479  if ( !fTCardCorrCellsNew[j] ) continue;
480 
481  // Accept only if at least 10 MeV
482  if ( fTCardCorrCellsEner[j] < 0.01 ) continue;
483 
484  // Add new cell
485  cellNumber = nCells;
486  absId = j;
487  amp = fTCardCorrCellsEner[j];
488  time = 615.*1e-9;
489  mclabel = -1;
490  efrac = 0.;
491 
492  Int_t ok = fCaloCells->SetCell(cellNumber, absId, amp, time, mclabel, efrac,1);
493 
494  if ( !ok ) AliError("Induced new cell could not be added!");
495 
496  nCells++;
497  }
498 
499  if ( nCellsNew > 0 ) fCaloCells->Sort();
500 }
501 
506 {
507  for(Int_t j = 0; j < fgkNEMCalCells; j++)
508  {
509  fTCardCorrCellsEner[j] = 0.;
510  fTCardCorrCellsNew [j] = kFALSE;
511  }
512 
513  if ( fAODCellsTmp ) fAODCellsTmp->DeleteContainer();
514 }
515 
525 {
526  if ( absID < 0 || absID >= 24*48*fGeom->GetNumberOfSuperModules() )
527  return kFALSE;
528 
529  Int_t imod = -1,iTower = -1, iIphi = -1, iIeta = -1;
530  if (!fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
531  return kFALSE;
532 
533  return kTRUE;
534 }
535 
540 {
541  printf("T-Card emulation activated, energy conservation <%d>, randomize E <%d>, induced energy parameters:\n",
543 
544  printf("T-Card emulation super-modules fraction: Min cell E %2.1f MeV; "
545  "induced Min E %2.1f MeV; Max at low E %2.1f MeV; Max E %2.2f GeV\n",
548 
549  for(Int_t ism = 0; ism < fgkNsm; ism++)
550  {
551  printf("\t sm %d, fraction %2.3f, E frac abs min %2.3e max %2.3e \n",
553 
554  for(Int_t icell = 0; icell < 4; icell++)
555  {
556  printf("\t \t cell type %d, c %2.4e, p0 %2.4e, p1 %2.4e, sigma %2.4e \n",
557  icell,fTCardCorrInduceEner[icell][ism],fTCardCorrInduceEnerFrac[icell][ism],
559  }
560  }
561 }
Float_t fTCardCorrInduceEnerFracP1[4][fgkNsm]
Induced energy loss gauss fraction param1 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
Bool_t fTCardCorrCellsNew[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells, that before had no signal.
AliEMCALGeometry * fGeom
! Geometry object
Float_t fTCardCorrMaxInduced
Maximum induced energy signal on adjacent cells.
double Double_t
Definition: External.C:58
static const Int_t fgkNsm
Total number of super-modules.
Float_t fTCardCorrMinInduced
Minimum induced energy signal on adjacent cells, sum of induced plus original energy, use same as cell energy clusterization cut.
static const Int_t fgkNEMCalCells
Total number of cells in the calorimeter, 10*48*24 (EMCal) + 4*48*8 (EMCal/DCal 1/3) + 6*32*24 (DCal)...
Correction component to emulate cell-level crosstalk in the EMCal correction framework.
Bool_t fTCardCorrClusEnerConserv
When making correlation, subtract from the reference cell the induced energy on the neighbour cells...
Float_t fTCardCorrInduceEnerProb[fgkNsm]
Probability to induce energy loss per SM.
void SetProperty(Float_t val[][fgkNsm], std::vector< double > &property, unsigned int iSM, const std::string &name)
AliAODCaloCells * fAODCellsTmp
! Temporal array of cells copy
AliVCaloCells * fCaloCells
! Pointer to CaloCells
Some utilities for cluster and cell treatment.
Float_t fTCardCorrMaxInducedELeak
Maximum value of induced energy signal that is always leaked, ~5-10 MeV.
AliEMCALRecoUtils * fRecoUtils
Pointer to RecoUtils.
Float_t fTCardCorrInduceEnerFrac[4][fgkNsm]
Induced energy loss gauss fraction param0 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
int Int_t
Definition: External.C:63
Float_t fTCardCorrMinAmp
Minimum cell energy to induce signal on adjacent cells.
Float_t fTCardCorrInduceEnerFracWidth[4][fgkNsm]
Induced energy loss gauss witdth on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ ...
float Float_t
Definition: External.C:68
Base class for correction components in the EMCal correction framework.
Float_t fTCardCorrCellsEner[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells.
Float_t fTCardCorrInduceEnerFracMin[fgkNsm]
In case fTCardCorrInduceEnerFracP1 is non null, restrict the minimum fraction of induced energy per S...
TList * fOutput
! List of output histograms
short Short_t
Definition: External.C:23
Bool_t fCreateHisto
Flag to make some basic histograms.
TH1F * fCellEnergyDistBefore
! cell energy distribution, before energy smearing
TH1F * fCellEnergyDistAfter
! cell energy distribution, after energy smearing
Float_t fTCardCorrInduceEner[4][fgkNsm]
Induced energy loss gauss constant on 0-same row, diff col, 1-up/down cells left/right col 2-left/rig...
static RegisterCorrectionComponent< AliEmcalCorrectionCellEmulateCrosstalk > reg
AliEmcalCorrectionEventManager fEventManager
Minimal task which inherits from AliAnalysisTaskSE and manages access to the event.
Float_t fTCardCorrInduceEnerFracMax[fgkNsm]
In case fTCardCorrInduceEnerFracP1 is non null, restrict the maximum fraction of induced energy per S...
bool Bool_t
Definition: External.C:53
void CalculateInducedEnergyInTCardCell(Int_t absId, Int_t absIdRef, Int_t sm, Float_t ampRef, Int_t cellCase)
bool GetProperty(std::string propertyName, T &property, bool requiredProperty=true, std::string correctionName="")
Retrieve property.