AliPhysics  fffcdf3 (fffcdf3)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEMCALClusterize.cxx
Go to the documentation of this file.
1  /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- Root ---
17 #include <TString.h>
18 #include <TRefArray.h>
19 #include <TClonesArray.h>
20 #include <TTree.h>
21 #include <TGeoManager.h>
22 #include <TROOT.h>
23 #include <TInterpreter.h>
24 #include <TFile.h>
25 
26 // --- AliRoot Analysis Steering
27 #include "AliAnalysisTask.h"
28 #include "AliAnalysisManager.h"
29 #include "AliESDEvent.h"
30 #include "AliGeomManager.h"
31 #include "AliVCaloCells.h"
32 #include "AliAODCaloCluster.h"
33 #include "AliCDBManager.h"
34 #include "AliCDBStorage.h"
35 #include "AliCDBEntry.h"
36 #include "AliLog.h"
37 #include "AliVEventHandler.h"
38 #include "AliAODInputHandler.h"
39 #include "AliOADBContainer.h"
40 #include "AliAODMCParticle.h"
41 #include "AliCentrality.h"
42 #include "AliMultSelection.h"
43 
44 // --- EMCAL
45 #include "AliEMCALAfterBurnerUF.h"
46 #include "AliEMCALGeometry.h"
47 #include "AliEMCALClusterizerNxN.h"
48 #include "AliEMCALClusterizerv1.h"
49 #include "AliEMCALClusterizerv2.h"
50 #include "AliEMCALRecPoint.h"
51 #include "AliEMCALDigit.h"
52 
54 
58 
59 //______________________________________________________________________________
60 // Constructor.
61 //______________________________________________________________________________
63 : AliAnalysisTaskSE(name)
64 , fEvent(0)
65 , fGeom(0), fGeomName("")
66 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
67 , fOCDBpath(""), fAccessOCDB(kFALSE)
68 , fDigitsArr(0), fClusterArr(0)
69 , fCaloClusterArr(0), fCaloCells(0)
70 , fRecParam(0), fClusterizer(0)
71 , fUnfolder(0), fJustUnfold(kFALSE)
72 , fOutputAODBranch(0), fOutputAODBranchName("")
73 , fOutputAODCells (0), fOutputAODCellsName ("")
74 , fOutputAODBranchSet(0)
75 , fFillAODFile(kFALSE), fFillAODHeader(0)
76 , fFillAODCaloCells(0), fRun(-1)
77 , fRecoUtils(0), fConfigName("")
78 , fOrgClusterCellId()
79 , fCellLabels(), fCellSecondLabels(), fCellTime()
80 , fCellMatchdEta(), fCellMatchdPhi()
81 , fRecalibrateWithClusterTime(0)
82 , fMaxEvent(0), fDoTrackMatching(kFALSE), fUpdateCell(0)
83 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
84 , fRejectBelowThreshold(kFALSE)
85 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
86 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
87 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
88 , fConstantTimeShift(0)
89 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
90 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
91 , fSetCellMCLabelFromCluster(0)
92 , fSetCellMCLabelFromEdepFrac(0)
93 , fRemapMCLabelForAODs(0)
94 , fInputFromFilter(0)
95 , fTCardCorrEmulation(0), fTCardCorrClusEnerConserv(0)
96 , fRandom(0), fRandomizeTCard(1)
97 , fTCardCorrMinAmp(0.01), fTCardCorrMaxInduced(100)
98 , fPrintOnce(0)
99 
100 {
101  for(Int_t i = 0; i < 22; i++)
102  {
103  fGeomMatrix[i] = 0;
105 
106  for(Int_t j = 0; j < 4 ; j++)
107  {
108  fTCardCorrInduceEner [j][i] = 0 ;
109  fTCardCorrInduceEnerFrac [j][i] = 0 ;
110  fTCardCorrInduceEnerFracP1 [j][i] = 0 ;
111  fTCardCorrInduceEnerFracWidth[j][i] = 0 ;
112  }
113  }
114 
115  ResetArrays();
116 
117  fCentralityBin[0] = fCentralityBin[1]=-1;
118 }
119 
120 //______________________________________________________________
122 //______________________________________________________________
124 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
125 , fEvent(0)
126 , fGeom(0), fGeomName("")
127 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
128 , fOCDBpath(""), fAccessOCDB(kFALSE)
129 , fDigitsArr(0), fClusterArr(0)
130 , fCaloClusterArr(0), fCaloCells(0)
131 , fRecParam(0), fClusterizer(0)
132 , fUnfolder(0), fJustUnfold(kFALSE)
133 , fOutputAODBranch(0), fOutputAODBranchName("")
134 , fOutputAODCells (0), fOutputAODCellsName ("")
135 , fOutputAODBranchSet(0)
136 , fFillAODFile(kFALSE), fFillAODHeader(0)
137 , fFillAODCaloCells(0), fRun(-1)
138 , fRecoUtils(0), fConfigName("")
139 , fOrgClusterCellId()
140 , fCellLabels(), fCellSecondLabels(), fCellTime()
141 , fCellMatchdEta(), fCellMatchdPhi()
142 , fRecalibrateWithClusterTime(0)
143 , fMaxEvent(0), fDoTrackMatching(kFALSE), fUpdateCell(0)
144 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
145 , fRejectBelowThreshold(kFALSE)
146 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
147 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
148 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
149 , fConstantTimeShift(0)
150 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
151 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
152 , fSetCellMCLabelFromCluster(0)
153 , fSetCellMCLabelFromEdepFrac(0)
154 , fRemapMCLabelForAODs(0)
155 , fInputFromFilter(0)
156 , fTCardCorrEmulation(0), fTCardCorrClusEnerConserv(0)
157 , fRandom(0), fRandomizeTCard(1)
158 , fTCardCorrMinAmp(0.01), fTCardCorrMaxInduced(100)
159 , fPrintOnce(0)
160 {
161  for(Int_t i = 0; i < 22; i++)
162  {
163  fGeomMatrix[i] = 0;
165 
166  for(Int_t j = 0; j < 4 ; j++)
167  {
168  fTCardCorrInduceEner [j][i] = 0 ;
169  fTCardCorrInduceEnerFrac [j][i] = 0 ;
170  fTCardCorrInduceEnerFracP1 [j][i] = 0 ;
171  fTCardCorrInduceEnerFracWidth[j][i] = 0 ;
172  }
173  }
174 
175  ResetArrays();
176 
177  fCentralityBin[0] = fCentralityBin[1]=-1;
178 }
179 
180 //_______________________________________________________________
182 //_______________________________________________________________
184 {
185  if (fDigitsArr)
186  {
187  fDigitsArr->Clear("C");
188  delete fDigitsArr;
189  }
190 
191  if (fClusterArr)
192  {
193  fClusterArr->Delete();
194  delete fClusterArr;
195  }
196 
197  if (fCaloClusterArr)
198  {
199  fCaloClusterArr->Delete();
200  delete fCaloClusterArr;
201  }
202 
203  if(fClusterizer) delete fClusterizer;
204  if(fUnfolder) delete fUnfolder;
205  if(fRecoUtils) delete fRecoUtils;
206 }
207 
208 //_______________________________________________________
216 //_______________________________________________________________________________
218 {
219  if ( absID < 0 || absID >= 24*48*fGeom->GetNumberOfSuperModules() )
220  return kFALSE;
221 
222  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
223  if (!fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
224  return kFALSE;
225 
226  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
227 
228  // Do not include bad channels found in analysis,
229  if ( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
230  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi) )
231  return kFALSE;
232 
233  return kTRUE;
234 }
235 
236 //_______________________________________________________
240 //_______________________________________________________
242 {
243  if(!fSelectEMCALEvent) return kTRUE; // accept
244 
245  if(fEMCALEnergyCut <= 0) return kTRUE; // accept
246 
247  Int_t nCluster = fEvent -> GetNumberOfCaloClusters();
248  Int_t bc = fEvent -> GetBunchCrossNumber();
249 
250  for(Int_t icalo = 0; icalo < nCluster; icalo++)
251  {
252  AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
253 
254  if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
255  fRecoUtils->IsGoodCluster(clus,fGeom,fCaloCells,bc))
256  {
257 
258  AliDebug(1, Form("Accept : E %2.2f > %2.2f, nCells %d > %d",
259  clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut));
260 
261  return kTRUE;
262  }
263 
264  }// loop
265 
266  AliDebug(1,"Reject");
267 
268  return kFALSE;
269 }
270 
271 //_______________________________________________
274 //_______________________________________________
276 {
277  // Set it only once, unless run changed
278  if ( fOADBSet ) return ;
279 
280  TString pass = GetPass();
281 
282  AliInfo(Form("Get AODB parameters from EMCAL in %s for run %d, and <%s>",fOADBFilePath.Data(),fRun,pass.Data()));
283 
284  Int_t nSM = fGeom->GetNumberOfSuperModules();
285 
286  // Bad map
287  if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
288  {
289  AliOADBContainer *contBC=new AliOADBContainer("");
290  contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
291 
292  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRun);
293 
294  if(arrayBC)
295  {
296  fRecoUtils->SwitchOnDistToBadChannelRecalculation();
297  AliInfo("Remove EMCAL bad cells");
298 
299  for (Int_t i=0; i < nSM; ++i)
300  {
301  TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
302 
303  if (hbm)
304  delete hbm;
305 
306  hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
307 
308  if (!hbm)
309  {
310  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
311  continue;
312  }
313 
314  hbm->SetDirectory(0);
315  fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
316 
317  } // loop
318  } else AliInfo("Do NOT remove EMCAL bad channels"); // run array
319 
320  delete contBC;
321  } // Remove bad
322 
323  // Energy Recalibration
324  if(fRecoUtils->IsRecalibrationOn())
325  {
326  AliOADBContainer *contRF=new AliOADBContainer("");
327 
328  contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
329 
330  TObjArray *recal=(TObjArray*)contRF->GetObject(fRun);
331 
332  if(recal)
333  {
334  TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
335 
336  if(recalpass)
337  {
338  TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
339 
340  if(recalib)
341  {
342  AliInfo("Recalibrate EMCAL");
343  for (Int_t i=0; i<nSM; ++i)
344  {
345  TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
346 
347  if (h)
348  delete h;
349 
350  h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
351 
352  if (!h)
353  {
354  AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
355  continue;
356  }
357 
358  h->SetDirectory(0);
359 
360  fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
361  } // SM loop
362  } else AliInfo("Do NOT recalibrate EMCAL, no params object array"); // array ok
363  } else AliInfo("Do NOT recalibrate EMCAL, no params for pass"); // array pass ok
364  } else AliInfo("Do NOT recalibrate EMCAL, no params for run"); // run number array ok
365 
366  delete contRF;
367  } // Recalibration on
368 
369  // Energy Recalibration, apply on top of previous calibration factors
370  if(fRecoUtils->IsRunDepRecalibrationOn())
371  {
372  AliOADBContainer *contRFTD=new AliOADBContainer("");
373 
374  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
375 
376  TH1S *htd=(TH1S*)contRFTD->GetObject(fRun);
377 
378  //If it did not exist for this run, get closes one
379  if (!htd)
380  {
381  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",fRun));
382  // let's get the closest fRun instead then..
383  Int_t lower = 0;
384  Int_t ic = 0;
385  Int_t maxEntry = contRFTD->GetNumberOfEntries();
386 
387  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < fRun) ) {
388  lower = ic;
389  ic++;
390  }
391 
392  Int_t closest = lower;
393  if ( (ic<maxEntry) &&
394  (contRFTD->LowerLimit(ic)-fRun) < (fRun - contRFTD->UpperLimit(lower)) ) {
395  closest = ic;
396  }
397 
398  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
399  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
400  }
401 
402  if(htd)
403  {
404  AliInfo("Recalibrate (Temperature) EMCAL");
405 
406  for (Int_t ism=0; ism<nSM; ++ism)
407  {
408  for (Int_t icol=0; icol<48; ++icol)
409  {
410  for (Int_t irow=0; irow<24; ++irow)
411  {
412  Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
413 
414  Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
415  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
416  //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
417  // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
418  fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
419  } // columns
420  } // rows
421  } // SM loop
422  } else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
423 
424  delete contRFTD;
425  } // Run by Run T calibration
426 
427  // Time Recalibration
428  if(fRecoUtils->IsTimeRecalibrationOn())
429  {
430  AliOADBContainer *contTRF=new AliOADBContainer("");
431 
432  contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
433 
434  TObjArray *trecal=(TObjArray*)contTRF->GetObject(fRun);
435 
436  if(trecal)
437  {
438  // pass number should be pass1 except on Run1 and special cases
439  TString passM = pass;
440  if ( pass=="spc_calo" ) passM = "pass3";
441  if ( fRun > 209121 ) passM = "pass1"; // run2 periods
442  if ( pass == "muon_calo_pass1" && fRun > 209121 && fRun < 244284 )
443  passM = "pass0";//period LHC15a-m
444 
445  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
446 
447  if(trecalpass)
448  {
449  AliInfo("Time Recalibrate EMCAL");
450  for (Int_t ibc = 0; ibc < 4; ++ibc)
451  {
452  TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
453 
454  if (h)
455  delete h;
456 
457  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
458 
459  if (!h)
460  {
461  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
462  continue;
463  }
464 
465  h->SetDirectory(0);
466 
467  fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
468  } // bunch crossing loop
469  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
470  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
471 
472  delete contTRF;
473  } // Time recalibration on
474 
475  // L1 Phase Time Recalibration
476  if( fRecoUtils->IsL1PhaseInTimeRecalibrationOn() )
477  {
478  // init default maps first
479  if (!fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationArray())
480  fRecoUtils->InitEMCALL1PhaseInTimeRecalibration() ;
481 
482  AliOADBContainer *contBC = new AliOADBContainer("");
483 
484  TFile *timeFile=new TFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"read");
485  if (!timeFile || timeFile->IsZombie())
486  {
487  AliFatal(Form("EMCALTimeL1PhaseCalib.root was not found in the path provided: %s",fOADBFilePath.Data()));
488  return ;
489  }
490 
491  if (timeFile) delete timeFile;
492 
493  contBC->InitFromFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"AliEMCALTimeL1PhaseCalib");
494 
495  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRun);
496  if (!arrayBC)
497  {
498  AliError(Form("No external L1 phase in time calibration set for run number: %d", fRun));
499  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
500  }
501  else
502  {
503  // Only 1 L1 phase correction possible, except special cases
504  TString pass2 = "pass1";
505 
506  if ( pass=="muon_calo_pass1" && fRun > 209121 && fRun < 244284 )
507  pass2 = "pass0"; // period LHC15a-m
508 
509  TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(pass2);
510  if (!arrayBCpass)
511  {
512  AliError(Form("No external L1 phase in time calibration set for: %d -%s", fRun,pass2.Data()));
513  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
514  }
515  else AliInfo("Recalibrate L1 Phase time");
516 
517  if(arrayBCpass)
518  {
519  if ( DebugLevel()>0 ) arrayBCpass->Print();
520 
521  TH1C *h = fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationForAllSM();
522  if (h) delete h;
523 
524  h = (TH1C*)arrayBCpass->FindObject(Form("h%d",fRun));
525 
526  if (!h)
527  {
528  AliFatal(Form("There is no calibration histogram h%d for this run",fRun));
529  return;
530  }
531 
532  h->SetDirectory(0);
533  fRecoUtils->SetEMCALL1PhaseInTimeRecalibrationForAllSM(h);
534  }
535  }
536 
537  delete contBC;
538  } // L1 Phase Time Recalibration
539 
540  // Parameters already set once, so do not it again, unless run changes
541  fOADBSet = kTRUE;
542 }
543 
544 //_________________________________________________
547 //_________________________________________________
549 {
550  // Set once per run
551  if ( fOADBSet ) return;
552 
553  AliDebug(1,"Begin");
554 
555  AliCDBManager *cdb = AliCDBManager::Instance();
556 
557  if (fOCDBpath.Length())
558  {
559  cdb->SetDefaultStorage(fOCDBpath.Data());
560  AliInfo(Form("Default storage %s",fOCDBpath.Data()));
561  }
562 
563  cdb->SetRun(fRun);
564 
565  //
566  // EMCAL from RAW OCDB
567  if (fOCDBpath.Contains("alien:"))
568  {
569  cdb->SetSpecificStorage("EMCAL/Calib/Data","raw://");
570  cdb->SetSpecificStorage("EMCAL/Calib/Time","raw://");
571  cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","raw://");
572  }
573 
574  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
575 
576  fOADBSet = kTRUE;
577 
578  return ;
579 }
580 
581 //_____________________________________________________
585 //_____________________________________________________
587 {
588  // Get the number of stored digits, to assign the index of the new ones
589  Int_t idigit = fDigitsArr->GetEntriesFast();
590 
591  Double_t fixTime = 615.*1e-9;
592 
593  for(Int_t j = 0; j < fgkNEMCalCells; j++)
594  {
595  // Newly created?
596  if ( !fTCardCorrCellsNew[j] ) continue;
597 
598  // Accept only if at least 10 MeV
599  if ( fTCardCorrCellsEner[j] < 0.01 ) continue;
600 
601  // Check if it was not masked
602  Float_t amp = 0;
603  Double_t time = 0;
604  if ( !fRecoUtils->AcceptCalibrateCell(j,0,amp,time,fCaloCells) ) continue;
605 
606 // printf("add new digit absId %d, accept? %d, digit %d, induced amp %2.2f\n",
607 // j,fRecoUtils->AcceptCalibrateCell(j,0,amp,time,fCaloCells),idigit,fTCardCorrCellsEner[j]);
608 
609  // Now add the cell to the digits list
610  //AliEMCALDigit* digit =
611  new((*fDigitsArr)[idigit]) AliEMCALDigit( -2, -2, j, fTCardCorrCellsEner[j], fixTime,AliEMCALDigit::kHG,idigit, 0, 0, 0);
612 
613  //digit->SetListOfParents(0,0x0,0x0); // not needed
614 
615  idigit++;
616  }// loop on all possible cells
617 }
618 
619 
620 //_____________________________________________________
625 //_____________________________________________________
627 {
628  fEvent = 0x0;
629 
630  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
631  Int_t eventN = Entry();
632  if(aodIH) eventN = aodIH->GetReadEntry();
633 
634  if (eventN > fMaxEvent)
635  return ;
636 
637  //printf("Clusterizer --- Event %d-- \n",eventN);
638 
639  // Check if input event are embedded events
640  // If so, take output event
641  if (aodIH && aodIH->GetMergeEvents())
642  {
643  fEvent = AODEvent();
644 
645  if(!aodIH->GetMergeEMCALCells())
646  AliFatal("Events merged but not EMCAL cells, check analysis settings!");
647 
648  AliDebug(1,"Use embedded events");
649 
650  AliDebug(1,Form("\t InputEvent N Clusters %d, N Cells %d",
651  InputEvent()->GetNumberOfCaloClusters(),InputEvent()->GetEMCALCells()->GetNumberOfCells()));
652 
653  AliDebug(1,Form("\t MergedEvent N Clusters %d, N Cells %d",
654  aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()));
655 
656 // if(DebugLevel() > 1)
657 // {
658 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
659 // {
660 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
661 // if(sigCluster->IsEMCAL()) AliInfo(Form("\t \t Signal cluster: i %d, E %f",icl,sigCluster->E()));
662 // }
663 // }
664 
665  AliDebug(1,Form("\t OutputEvent N Clusters %d, N Cells %d",
666  AODEvent()->GetNumberOfCaloClusters(), AODEvent()->GetEMCALCells()->GetNumberOfCells()));
667  }
668  else if(fInputFromFilter)
669  {
670  //printf("Get Input From Filtered AOD\n");
671  fEvent = AODEvent();
672  }
673  else
674  {
675  fEvent = InputEvent();
678  }
679 
680  if (!fEvent)
681  {
682  AliError("Event not available");
683  return ;
684  }
685 
686  //Recover the pointer to CaloCells container
687  fCaloCells = fEvent->GetEMCALCells();
688 
689  //Process events if there is a high energy cluster
690  if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
691 
692  //-------------------------------------------------------------------------------------
693  // Reject events if LED was firing, use only for LHC11a data
694  // Reject event if triggered by exotic cell and remove exotic cells if not triggered
695  //-------------------------------------------------------------------------------------
696 
697  if( IsLEDEvent( fRun ) ) { fEvent = 0x0 ; return ; }
698 
699  if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
700 
701  //-------------------------------------------------------------------------------------
702  // Set the cluster array in the event (output or input)
703  //-------------------------------------------------------------------------------------
704 
705  if ( fFillAODFile )
706  {
707  //Magic line to write events to AOD filem put after event rejection
708  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
709  }
710  else if( !fOutputAODBranchSet )
711  {
712  // Create array of clusters/cells and put it in the input event, if output AOD not selected, only once
713  InputEvent()->AddObject(fOutputAODBranch);
714  AliInfo(Form("Add AOD clusters branch <%s> to input event",fOutputAODBranchName.Data()));
715 
716  if ( fOutputAODBranchName.Length() > 0 )
717  {
718  InputEvent()->AddObject(fOutputAODCells);
719  AliInfo(Form("Add AOD cells branch <%s> to input event",fOutputAODCellsName.Data()));
720  }
721 
722  fOutputAODBranchSet = kTRUE;
723  }
724 }
725 
726 //____________________________________________________
731 //____________________________________________________
733 {
734  Int_t nClusters = fEvent->GetNumberOfCaloClusters();
735  Int_t nClustersOrg = 0;
736 
737  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
738  if(aodIH && aodIH->GetEventToMerge()) //Embedding
739  nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
740 
741  ResetArrays();
742 
743  // Loop on original clusters, get MC labels, cluster time (OLD AODs),
744  // or track matching residuals (if matching is not requested)
746  {
747  for (Int_t i = 0; i < nClusters; i++)
748  {
749  AliVCluster *clus = 0;
750  if(aodIH && aodIH->GetEventToMerge()) //Embedding
751  clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
752  else
753  clus = fEvent->GetCaloCluster(i);
754 
755  if(!clus) return;
756 
757  nClustersOrg++;
758 
759  if(!clus->IsEMCAL()) continue;
760 
761  Int_t label = clus->GetLabel();
762  Int_t label2 = -1 ;
763  if (clus->GetNLabels() >=2 ) label2 = clus->GetLabelAt(1) ;
764 
765  //printf("Org cluster %d) ID = %d, E %2.2f, Time %3.0f, N Cells %d, N MC labels %d, main MC label %d, all MC labels:\n",
766  // i, clus->GetID(), clus->E(), clus->GetTOF()*1e9, clus->GetNCells(), clus->GetNLabels(), label);
767  //
768  //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++)
769  // printf("%d) Label %d, E dep frac %0.2f; ",
770  // imc, clus->GetLabelAt(imc),clus->GetClusterMCEdepFraction(imc));
771  //if(clus->GetNLabels() > 0) printf("\n");
772 
773  UShort_t * index = clus->GetCellsAbsId() ;
774  for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
775  {
776  fOrgClusterCellId[index[icell]] = i;
777  fCellTime[index[icell]] = clus->GetTOF();
778  fCellMatchdEta[index[icell]] = clus->GetTrackDz();
779  fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
780 
782  {
783  fCellLabels[index[icell]] = label;
784  fCellSecondLabels[index[icell]] = label2;
785  }
786  }
787  }
788  }
789 
790  // Do here induced cell energy assignation by T-Card correlation emulation, ONLY MC
792 
793  //
794  // Transform CaloCells into Digits
795  //
796  Int_t idigit = 0;
797  Int_t id = -1;
798  Float_t amp = -1;
799  Double_t time = -1;
800 
801  Int_t bc = InputEvent()->GetBunchCrossNumber();
802 
803  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
804  {
805  // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
806  id = fCaloCells->GetCellNumber(icell);
807  Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,fCaloCells);
808  time-=fConstantTimeShift*1e-9; // only in case of simulations done before 2015
809 
810  // Do not include cells with too low energy, nor exotic cell
811  if( amp < fRecParam->GetMinECut() ||
812  time > fRecParam->GetTimeMax() ||
813  time < fRecParam->GetTimeMin() ) accept = kFALSE;
814 
815  // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
817  {
818  time = fCellTime[id];
819  //printf("cell %d time org %f - ",id, time*1.e9);
820  fRecoUtils->RecalibrateCellTime(id,bc,time);
821  //printf("recal %f\n",time*1.e9);
822  }
823 
824  //Exotic?
825  if (accept && fRecoUtils->IsExoticCell(id,fCaloCells,bc))
826  accept = kFALSE;
827 
828  if( !accept )
829  {
830  AliDebug(2,Form("Remove channel absId %d, index %d of %d, amp %f, time %f",
831  id,icell, fCaloCells->GetNumberOfCells(), amp, time*1.e9));
832  continue;
833  }
834 
835  //
836  // MC
837  //
838 
839  // Old way
840  Int_t mcLabel = fCaloCells->GetMCLabel(icell);
841  Float_t eDep = amp;
842 
843  //printf("--- Selected cell %d--- amp %f\n",id,amp);
844 
845  // New way
846  TArrayI labeArr(0);
847  TArrayF eDepArr(0);
848  Int_t nLabels = 0;
849 
851  {
852  // Old way to recover/set the cell MC label
853  // Only possibility for old Run1 productions
854 
855  if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
856 
857  else if( fSetCellMCLabelFromCluster == 0 &&
859 
860  else mcLabel = -1; // found later
861 
862  // Last parameter of the digit object should be MC deposited energy,
863  // since it is not available in aliroot before year 2016, add just the cell amplitude so that
864  // we give more weight to the MC label of the cell with highest energy in the cluster
865 
866  Float_t efrac = fCaloCells->GetEFraction(icell);
867 
868  // When checking the MC of digits, give weight to cells with embedded signal
869  if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
870 
871  eDep *= efrac ;
872  }
873  else // fSetCellMCLabelFromEdepFrac = true
874  {
875  // New way, valid only for MC productions with aliroot > v5-07-21
876  mcLabel = -1;
877 
878  // Map the digit to cell index for later to calculate the cell MC energy deposition map
879  fCellLabels[id] = idigit;
880  //printf("\t absId %d, idigit %d\n",id,idigit);
881 
882  if(fOrgClusterCellId[id] < 0) continue; // index can be negative if noisy cell that did not form cluster
883 
884  AliVCluster *clus = 0;
885  Int_t iclus = fOrgClusterCellId[id];
886 
887  if(iclus < 0)
888  {
889  AliInfo("Negative original cluster index, skip \n");
890  continue;
891  }
892 
893  if(aodIH && aodIH->GetEventToMerge()) //Embedding
894  clus = aodIH->GetEventToMerge()->GetCaloCluster(iclus); //Get clusters directly from embedded signal
895  else
896  clus = fEvent->GetCaloCluster(iclus);
897 
898  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(id, clus, MCEvent(), amp, labeArr, eDepArr);
899  nLabels = labeArr.GetSize();
900 
901  } // cell MC label, new
902 
903  // Apply here the found induced energies
904  if ( fTCardCorrEmulation )
905  {
906 // if( TMath::Abs(fTCardCorrCellsEner[id]) > 0.001 )
907 // printf("add energy to digit %d, absId %d: amp %2.2f + %2.2f\n",idigit,id,amp,fTCardCorrCellsEner[id]);
908  amp+=fTCardCorrCellsEner[id];
909  }
910 
911  //
912  // Create the digit
913  //
914  if(amp <= 0.01) continue ; // accept if > 10 MeV
915 
916  //if(amp > 1.5)printf("*** Add digit *** amp %f, nlabels %d, label %d, found %d, edep fract tot %f, ncluslabels %d\n",amp, nLabels,mcLabel,found,edepTotFrac,nClusLabels);
917 
918  AliEMCALDigit* digit = new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, eDep);
919 
920  if(nLabels > 0)
921  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
922 
923  idigit++;
924 
925  }
926 
927  fDigitsArr->Sort();
928 
929  // Call after Sort, if not it screws up digits index order in clusterization
931 
932 
933  //-------------------------------------------------------------------------------------
934  // Do the clusterization
935  //-------------------------------------------------------------------------------------
936 
937  fClusterizer->Digits2Clusters("");
938 
939  //-------------------------------------------------------------------------------------
940  // Transform the recpoints into AliVClusters
941  //-------------------------------------------------------------------------------------
942 
944 
945  if(!fCaloClusterArr)
946  {
947  AliWarning(Form("No array with CaloClusters, input RecPoints entries %d",fClusterArr->GetEntriesFast()));
948  return;
949  }
950 
951  AliDebug(1,Form("N clusters: before recluster %d, after recluster %d, recpoints %d",
952  nClustersOrg, fCaloClusterArr->GetEntriesFast(),fClusterArr->GetEntriesFast()));
953 
954 // if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
955 // {
956 // AliInfo("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
957 // fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
958 // fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
959 // }
960 }
961 
962 
963 //_____________________________________________________
965 //_____________________________________________________
967 {
968  Double_t cellAmplitude = 0;
969  Double_t cellTime = 0;
970  Short_t cellNumber = 0;
971  Int_t cellMCLabel = 0;
972  Double_t cellEFrac = 0;
973  Int_t nClustersOrg = 0;
974 
975  // Fill the array with the EMCAL clusters, copy them
976  for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
977  {
978  AliVCluster *clus = fEvent->GetCaloCluster(i);
979  if(clus->IsEMCAL())
980  {
981  //recalibrate/remove bad channels/etc if requested
982  if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
983  {
984  continue;
985  }
986 
987  if(fRecoUtils->IsRecalibrationOn())
988  {
989  //Calibrate cluster
990  fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, fCaloCells);
991 
992  //CalibrateCells
993  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
994  {
995  if (fCaloCells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
996  break;
997 
998  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
999  fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
1000  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1001 
1002  //Do not include bad channels found in analysis?
1003  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
1004  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
1005  continue;
1006 
1007  fCaloCells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
1008 
1009  }// cells loop
1010  }// recalibrate
1011 
1012  //Cast to ESD or AOD, needed to create the cluster array
1013  AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
1014  AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
1015 
1016  if (esdCluster)
1017  {
1018  fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
1019  }//ESD
1020  else if(aodCluster)
1021  {
1022  fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
1023  }//AOD
1024  else
1025  AliWarning("Wrong CaloCluster type?");
1026 
1027  nClustersOrg++;
1028  }
1029  }
1030 
1031  //Do the unfolding
1032  fUnfolder->UnfoldClusters(fCaloClusterArr, fCaloCells);
1033 
1034  //CLEAN-UP
1035  fUnfolder->Clear();
1036 }
1037 
1038 //_______________________________________________________________
1051 //_______________________________________________________________
1053 (Bool_t bMC , Bool_t bExotic, Bool_t bNonLin,
1054  Bool_t bRecalE, Bool_t bBad , Bool_t bRecalT, Int_t debug)
1055 {
1056  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - **** Start ***\n");
1057 
1058  // Init
1059  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils ;
1060 
1061  // Exotic cells removal
1062 
1063  if(bExotic)
1064  {
1065  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Remove exotics in EMCAL\n");
1066  fRecoUtils->SwitchOnRejectExoticCell() ;
1067  fRecoUtils->SwitchOnRejectExoticCluster();
1068 
1069  // fRecoUtils->SetExoticCellDiffTimeCut(50); // If |t cell max - t cell in cross| > 50 do not add its energy, avoid
1070  fRecoUtils->SetExoticCellFractionCut(0.97); // 1-Ecross/Ecell > 0.97 -> out
1071  fRecoUtils->SetExoticCellMinAmplitudeCut(4.); // 4 GeV
1072  }
1073 
1074  // Recalibration factors
1075 
1076  if(bRecalE && ! bMC)
1077  {
1078  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on energy recalibration in EMCAL\n");
1079  fRecoUtils->SwitchOnRecalibration();
1080  fRecoUtils->SwitchOnRunDepCorrection();
1081  }
1082 
1083  // Remove EMCAL hot channels
1084 
1085  if(bBad)
1086  {
1087  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on bad channels removal in EMCAL\n");
1088  fRecoUtils->SwitchOnBadChannelsRemoval();
1089  fRecoUtils->SwitchOnDistToBadChannelRecalculation();
1090  }
1091 
1092  // *** Time recalibration settings ***
1093 
1094  if(bRecalT && ! bMC)
1095  {
1096  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on time recalibration in EMCAL\n");
1097  fRecoUtils->SwitchOnTimeRecalibration();
1098  fRecoUtils->SwitchOnL1PhaseInTimeRecalibration() ;
1099  }
1100 
1101  // Recalculate position with method
1102 
1103  fRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
1104 
1105  // Non linearity
1106 
1107  if( bNonLin )
1108  {
1109  if(!bMC)
1110  {
1111  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kBeamTestCorrected xxx\n");
1112  fRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrectedv3);
1113  }
1114  else
1115  {
1116  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kPi0MCv3 xxx\n");
1117  fRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MCv3);
1118  }
1119  }
1120  else
1121  {
1122  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx DON'T SET Non linearity correction xxx\n");
1123  fRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kNoCorrection);
1124  }
1125 
1126 }
1127 
1128 //_____________________________________________________
1130 //_____________________________________________________
1132 {
1133  AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
1134  Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
1135 
1136  AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
1137  aodEMcells.CreateContainer(nEMcell);
1138  aodEMcells.SetType(AliVCaloCells::kEMCALCell);
1139  Double_t calibFactor = 1.;
1140  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1141  {
1142  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1143  fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
1144  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1145 
1146  if(fRecoUtils->IsRecalibrationOn())
1147  {
1148  calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
1149  }
1150 
1151  if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
1152  { //Channel is not declared as bad
1153  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
1154  eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
1155  }
1156  else
1157  {
1158  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
1159  }
1160  }
1161 
1162  aodEMcells.Sort();
1163 }
1164 
1165 //__________________________________________________
1167 //__________________________________________________
1169 {
1170  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
1171  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
1172 
1173  Double_t pos[3] ;
1174  Double_t covVtx[6];
1175  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
1176 
1177  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
1178  if(!header) AliFatal("Not a standard AOD");
1179  header->SetRunNumber(fRun);
1180 
1181  if(esdevent)
1182  {
1183  TTree* tree = fInputHandler->GetTree();
1184  if (tree)
1185  {
1186  TFile* file = tree->GetCurrentFile();
1187  if (file) header->SetESDFileName(file->GetName());
1188  }
1189  }
1190  else if (aodevent) {
1191  AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
1192  if(!aodheader) AliFatal("Not a standard AOD");
1193  header->SetESDFileName(aodheader->GetESDFileName());
1194  }
1195 
1196  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
1197  header->SetOrbitNumber(fEvent->GetOrbitNumber());
1198  header->SetPeriodNumber(fEvent->GetPeriodNumber());
1199  header->SetEventType(fEvent->GetEventType());
1200 
1201  // Centrality
1202  if(fUseAliCentrality)
1203  {
1204  if(fEvent->GetCentrality())
1205  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
1206  else
1207  header->SetCentrality(0);
1208  }
1209 
1210  // Trigger
1211  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
1212 
1213  header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
1214 
1215  header->SetTriggerMask(fEvent->GetTriggerMask());
1216  header->SetTriggerCluster(fEvent->GetTriggerCluster());
1217 
1218  if (esdevent)
1219  {
1220  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
1221  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
1222  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
1223  }
1224  else if (aodevent)
1225  {
1226  header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
1227  header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
1228  header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
1229  }
1230 
1231  header->SetMagneticField(fEvent->GetMagneticField());
1232  //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
1233 
1234  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
1235  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
1236  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
1237  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
1238  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
1239 
1240  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
1241  Float_t diamcov[3];
1242  fEvent->GetDiamondCovXY(diamcov);
1243  header->SetDiamond(diamxy,diamcov);
1244  if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
1245  else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
1246 
1247  //
1248  Int_t nVertices = 1 ;/* = prim. vtx*/;
1249  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1250 
1251  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1252 
1253  // Access to the AOD container of vertices
1254  TClonesArray &vertices = *(AODEvent()->GetVertices());
1255  Int_t jVertices=0;
1256 
1257  // Add primary vertex. The primary tracks will be defined
1258  // after the loops on the composite objects (V0, cascades, kinks)
1259  fEvent->GetPrimaryVertex()->GetXYZ(pos);
1260  Float_t chi = 0;
1261  if (esdevent)
1262  {
1263  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1264  chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
1265  }
1266  else if (aodevent)
1267  {
1268  aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1269  chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
1270  }
1271 
1272  AliAODVertex * primary = new(vertices[jVertices++])
1273  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
1274  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
1275  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
1276 }
1277 
1278 //___________________________________________________________
1282 //___________________________________________________________
1284 {
1285  // First recalculate track-matching for the new clusters
1286  if(fDoTrackMatching)
1287  {
1288  fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom,MCEvent());
1289  }
1290 
1291  // Put the new clusters in the AOD list
1292 
1293  Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
1294 
1295  for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
1296  {
1297  AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1298 
1299  newCluster->SetID(i);
1300 
1301  // Correct cluster energy non linearity
1302  newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
1303 
1304  //Add matched track
1305  if(fDoTrackMatching)
1306  {
1307  Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1308  if(trackIndex >= 0)
1309  {
1310  newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1311  AliDebug(2,Form("Matched Track index %d to new cluster %d",trackIndex,i));
1312  }
1313 
1314  Float_t dR = 999., dZ = 999.;
1315  fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dZ,dR);
1316  newCluster->SetTrackDistance(dR,dZ);
1317  }
1318  else
1319  {// Assign previously assigned matched track in reco, very very rough
1320  Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1321  newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1322  }
1323 
1324  //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1325  //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1326  //printf("\n");
1327  //
1328  //printf("New cluster %d) ID = %d, E %2.2f, Time %3.0f, N Cells %d, N MC labels %d, main MC label %d, all MC labels:\n",
1329  // i, newCluster->GetID(), newCluster->E(), newCluster->GetTOF()*1e9, newCluster->GetNCells(), newCluster->GetNLabels(), newCluster->GetLabel());
1330  //
1331  //for(Int_t imc = 0; imc < newCluster->GetNLabels(); imc++)
1332  // printf("%d) Label %d, E dep frac %0.2f; ",
1333  // imc,newCluster->GetLabelAt(imc),newCluster->GetClusterMCEdepFraction(imc));
1334  //if(newCluster->GetNLabels() > 0) printf("\n");
1335 
1336  // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1337  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, newCluster);
1338 
1339  new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
1340 
1341  AliDebug(2,Form("New cluster %d of %d, energy %f, mc label %d",
1342  newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()));
1343 
1344  } // cluster loop
1345 
1346  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1347 }
1348 
1349 
1350 //________________________________________________________________
1352 //________________________________________________________________
1354 {
1355  if(fUseAliCentrality)
1356  {
1357  if(GetCentrality()) return GetCentrality()->GetCentralityPercentile(fCentralityClass) ;
1358  else return -1. ;
1359  }
1360  else
1361  {
1362  if(GetMultSelCen()) return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass,kTRUE) ;
1363  else return -1. ;
1364  }
1365 
1366  return -1.;
1367 }
1368 
1369 //_______________________________________________
1371 //_______________________________________________
1373 {
1374  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1375  {
1376  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null");
1377  return TString("");
1378  }
1379 
1380  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1381  {
1382  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null");
1383  return TString("");
1384  }
1385 
1386  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1387  if (pass.Contains("ass1")) return TString("pass1");
1388  else if (pass.Contains("ass2")) return TString("pass2");
1389  else if (pass.Contains("ass3")) return TString("pass3");
1390  else if (pass.Contains("ass4")) return TString("pass4");
1391  else if (pass.Contains("ass5")) return TString("pass5");
1392  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1393  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1394  {
1395  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1396  return TString("pass1");
1397  }
1398  else if (pass.Contains("LHC14a1a"))
1399  {
1400  AliInfo("Enable EMCal energy calibration for this MC production!!");
1401 
1402  fRecoUtils->SwitchOnRecalibration();
1403 
1404  return TString("LHC14a1a");
1405  }
1406 
1407  // No condition fullfilled, give a default value
1408  AliInfo("Pass number string not found");
1409 
1410  return TString("");
1411 }
1412 
1413 //_________________________________________
1416 //_________________________________________
1418 {
1419  if(fDebug >=0) (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1420 
1421  fOADBSet = kFALSE;
1422  if(fOADBFilePath == "") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
1423 
1424  fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1425 
1426  if(!fRecParam) fRecParam = new AliEMCALRecParam;
1427  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1428 
1429  if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1430  if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1431  if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1432  fRejectBelowThreshold = kFALSE;
1433 
1434  // Centrality
1435  if(fCentralityClass == "") fCentralityClass = "V0M";
1436 
1437  if (fOCDBpath == "") fOCDBpath = "raw://" ;
1438  if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1439 
1440  if(gROOT->LoadMacro(fConfigName) >=0)
1441  {
1442  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
1443  AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1444  fGeomName = clus->fGeomName;
1446  fOCDBpath = clus->fOCDBpath;
1447  fAccessOCDB = clus->fAccessOCDB;
1448  fRecParam = clus->fRecParam;
1449  fJustUnfold = clus->fJustUnfold;
1450  fFillAODFile = clus->fFillAODFile;
1451  fRecoUtils = clus->fRecoUtils;
1452  fConfigName = clus->fConfigName;
1453  fMaxEvent = clus->fMaxEvent;
1455  fUpdateCell = clus->fUpdateCell;
1457  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1459  fCentralityBin[0] = clus->fCentralityBin[0];
1460  fCentralityBin[1] = clus->fCentralityBin[1];
1462  }
1463 }
1464 
1465 //_______________________________________________________
1468 //_______________________________________________________
1470 {
1471  if (fJustUnfold)
1472  {
1473  // init the unfolding afterburner
1474  delete fUnfolder;
1475  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1476  return;
1477  }
1478 
1479  //First init the clusterizer
1480  delete fClusterizer;
1481  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1482  fClusterizer = new AliEMCALClusterizerv1 (fGeom);
1483  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1484  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1485  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1486  {
1487  fClusterizer = new AliEMCALClusterizerNxN(fGeom);
1488  fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1489  fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1490  }
1491  else
1492  {
1493  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1494  }
1495 
1496  //Now set the parameters
1497  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1498  fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1499  fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1500  fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1501  fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1502  fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1503  fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1504  fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1505  fClusterizer->SetTimeCalibration ( fRecParam->IsTimeCalibrationOn() );
1506  fClusterizer->SetInputCalibrated ( kTRUE );
1507  fClusterizer->SetJustClusters ( kTRUE );
1508 
1509  // Initialize the cluster rec points and digits arrays and get them.
1510  fClusterizer->SetOutput(0);
1511  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1512  fDigitsArr = fClusterizer->GetDigits();
1513 
1514  //In case of unfolding after clusterization is requested, set the corresponding parameters
1515  if(fRecParam->GetUnfold())
1516  {
1517  Int_t i=0;
1518  for (i = 0; i < 8; i++)
1519  {
1520  fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1521  }//end of loop over parameters
1522 
1523  for (i = 0; i < 3; i++)
1524  {
1525  fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1526  fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1527  }//end of loop over parameters
1528 
1529  fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1530  fClusterizer->InitClusterUnfolding();
1531 
1532  }// to unfold
1533 }
1534 
1535 //________________________________________________________________
1539 //________________________________________________________________
1541 {
1542  if(fGeomMatrixSet) return;
1543 
1544  if (!fGeom)
1545  {
1546  if(fGeomName=="")
1547  {
1548  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(fRun);
1549  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fGeom->GetName(),fRun));
1550  }
1551  else
1552  {
1553  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1554  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeomName.Data()));
1555  }
1556 
1557  // Init geometry, I do not like much to do it like this ...
1558  if(fImportGeometryFromFile && !gGeoManager)
1559  {
1560  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1561  {
1562  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1563  if (fRun < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1564  else if (fRun < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1565  else if (fRun < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1566  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >=2015
1567  }
1568 
1569  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1570 
1571  TGeoManager::Import(fImportGeometryFilePath) ;
1572  }
1573 
1574  AliDebug(1,Form("Init for run=%d",fRun));
1575  if (!gGeoManager) AliDebug(1,"Careful!, gGeoManager not loaded, load misalign matrices");
1576  } // geometry pointer did not exist before
1577 
1578  if(fLoadGeomMatrices)
1579  {
1580  AliInfo("Load user defined EMCAL geometry matrices");
1581 
1582  // OADB if available
1583  AliOADBContainer emcGeoMat("AliEMCALgeo");
1584  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1585  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(fRun,"EmcalMatrices");
1586 
1587  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1588  {
1589 
1590  if (!fGeomMatrix[mod]) // Get it from OADB
1591  {
1592  AliDebug(2,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
1593  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1594 
1595  fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1596  }
1597 
1598  if(fGeomMatrix[mod])
1599  {
1600  if(DebugLevel() > 1)
1601  fGeomMatrix[mod]->Print();
1602 
1603  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1604  }
1605  else if(gGeoManager)
1606  {
1607  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
1608  fGeom->SetMisalMatrix(fGeom->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
1609  }
1610  else
1611  {
1612  AliError(Form("Alignment atrix for SM %d is not available",mod));
1613  }
1614 
1615  fGeomMatrixSet=kTRUE;
1616 
1617  }//SM loop
1618  }//Load matrices
1619  else if(!gGeoManager)
1620  {
1621  AliInfo("Get geo matrices from data");
1622  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1623  if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1624  {
1625  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
1626  }//AOD
1627  else
1628  {
1629  AliDebug(1,"Load Misaligned matrices");
1630 
1631  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1632 
1633  if(!esd)
1634  {
1635  AliError("This event does not contain ESDs?");
1636  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1637  return;
1638  }
1639 
1640  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1641  {
1642  if(DebugLevel() > 1)
1643  esd->GetEMCALMatrix(mod)->Print();
1644 
1645  if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1646 
1647  }
1648 
1649  fGeomMatrixSet=kTRUE;
1650 
1651  } // ESD
1652  } // Load matrices from Data
1653 }
1654 
1655 //____________________________________________________
1660 //____________________________________________________
1662 {
1663  if(!fRemoveExoticEvents) return kFALSE;
1664 
1665  // Loop on cells
1666 
1667  Float_t totCellE = 0;
1668  Int_t bc = InputEvent()->GetBunchCrossNumber();
1669 
1670  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1671  {
1672  Float_t ecell = 0 ;
1673  Double_t tcell = 0 ;
1674 
1675  Int_t absID = fCaloCells->GetCellNumber(icell);
1676  Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,fCaloCells);
1677  tcell-=fConstantTimeShift*1e-9;// Only for MC simulations done before 2015
1678 
1679  if(accept && !fRecoUtils->IsExoticCell(absID,fCaloCells,bc)) totCellE += ecell;
1680  }
1681 
1682  // TString triggerclasses = event->GetFiredTriggerClasses();
1683  // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1684  // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1685  // return;
1686  //
1687 
1688  //printf("TotE cell %f\n",totCellE);
1689  if(totCellE < 1) return kTRUE;
1690 
1691  return kFALSE;
1692 }
1693 
1694 //________________________________________________________________
1696 //________________________________________________________________
1698 {
1699  if(!fRemoveLEDEvents) return kFALSE;
1700 
1701  //check events of LHC11a period
1702  if(run < 146858 || run > 146860) return kFALSE ;
1703 
1704  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1705  Int_t ncellsSM3 = 0;
1706  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1707  {
1708  if ( fCaloCells->GetAmplitude (icell) > 0.1 &&
1709  fCaloCells->GetCellNumber(icell)/(24*48)==3 ) ncellsSM3++;
1710  }
1711 
1712  TString triggerclasses = fEvent->GetFiredTriggerClasses();
1713 
1714  Int_t ncellcut = 21;
1715  if(triggerclasses.Contains("EMC")) ncellcut = 35;
1716 
1717  if( ncellsSM3 >= ncellcut)
1718  {
1719  AliInfo(Form("Reject event %d with ncells in SM3 %d",(Int_t)Entry(),ncellsSM3));
1720  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1721  return kTRUE;
1722  }
1723 
1724  return kFALSE;
1725 }
1726 
1727 
1728 //_______________________________________________________
1731 //_______________________________________________________
1733 {
1734  Int_t id = -1;
1735  Float_t amp = -1;
1736 
1737  // Loop on all cells with signal
1738  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1739  {
1740  id = fCaloCells->GetCellNumber(icell);
1741  amp = fCaloCells->GetAmplitude (icell); // fCaloCells->GetCellAmplitude(id);
1742 
1743  if ( amp <= fTCardCorrMinAmp ) continue ;
1744 
1745  //
1746  // First get the SM, col-row of this tower
1747  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1748  fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
1749  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1750 
1751  //
1752  // Determine randomly if we want to create a correlation for this cell,
1753  // depending the SM number of the cell
1754  Float_t rand = fRandom.Uniform(0, 1);
1755 
1756  if ( rand > fTCardCorrInduceEnerProb[imod] ) continue;
1757 
1758  AliDebug(1,Form("Reference cell absId %d, iEta %d, iPhi %d, amp %2.3f",id,ieta,iphi,amp));
1759 
1760  //
1761  // Get the absId of the cells in the cross and same T-Card
1762  Int_t absIDup = -1;
1763  Int_t absIDdo = -1;
1764  Int_t absIDlr = -1;
1765  Int_t absIDuplr = -1;
1766  Int_t absIDdolr = -1;
1767 
1768  Int_t absIDup2 = -1;
1769  Int_t absIDup2lr = -1;
1770  Int_t absIDdo2 = -1;
1771  Int_t absIDdo2lr = -1;
1772 
1773  // Only 2 columns in the T-Card, +1 for even and -1 for odd with respect reference cell
1774  Int_t colShift = 0;
1775  if ( (ieta%2) && ieta <= AliEMCALGeoParams::fgkEMCALCols-1 ) colShift = -1;
1776  if ( !(ieta%2) && ieta >= 0 ) colShift = +1;
1777 
1778  absIDlr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+colShift);
1779 
1780  // Check up / down cells from reference cell not out of SM and in same T-Card
1781  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1 )
1782  {
1783  absIDup = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
1784  absIDuplr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta+colShift);
1785  }
1786 
1787  if ( iphi > 0 )
1788  {
1789  absIDdo = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
1790  absIDdolr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta+colShift);
1791  }
1792 
1793  // Check 2 up / 2 down cells from reference cell not out of SM
1794  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-2 )
1795  {
1796  absIDup2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta);
1797  absIDup2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta+colShift);
1798  }
1799 
1800  if ( iphi > 1 )
1801  {
1802  absIDdo2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta);
1803  absIDdo2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta+colShift);
1804  }
1805 
1806  // In same T-Card?
1807  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+1)/8) ) { absIDup = -1 ; absIDuplr = -1 ; }
1808  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-1)/8) ) { absIDdo = -1 ; absIDdolr = -1 ; }
1809  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+2)/8) ) { absIDup2 = -1 ; absIDup2lr = -1 ; }
1810  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-2)/8) ) { absIDdo2 = -1 ; absIDdo2lr = -1 ; }
1811 
1812  //
1813  // Check if they are not declared bad or exist
1814  Bool_t okup = AcceptCell(absIDup );
1815  Bool_t okdo = AcceptCell(absIDdo );
1816  Bool_t oklr = AcceptCell(absIDlr );
1817  Bool_t okuplr = AcceptCell(absIDuplr );
1818  Bool_t okdolr = AcceptCell(absIDdolr );
1819  Bool_t okup2 = AcceptCell(absIDup2 );
1820  Bool_t okdo2 = AcceptCell(absIDdo2 );
1821  Bool_t okup2lr= AcceptCell(absIDup2lr);
1822  Bool_t okdo2lr= AcceptCell(absIDdo2lr);
1823 
1824  AliDebug(1,Form("Same T-Card cells:\n \t up %d (%d), down %d (%d), left-right %d (%d), up-lr %d (%d), down-lr %d (%d)\n"
1825  "\t up2 %d (%d), down2 %d (%d), up2-lr %d (%d), down2-lr %d (%d)",
1826  absIDup ,okup ,absIDdo ,okdo ,absIDlr,oklr,absIDuplr ,okuplr ,absIDdolr ,okdolr ,
1827  absIDup2,okup2,absIDdo2,okdo2, absIDup2lr,okup2lr,absIDdo2lr,okdo2lr));
1828 
1829  //
1830  // Generate some energy for the nearby cells in same TCard , depending on this cell energy
1831  // Check if originally the tower had no or little energy, in which case tag it as new
1832  Float_t fracupdown = fTCardCorrInduceEnerFrac[0][imod]+amp*fTCardCorrInduceEnerFracP1[0][imod];
1833  Float_t fracupdownleri = fTCardCorrInduceEnerFrac[1][imod]+amp*fTCardCorrInduceEnerFracP1[1][imod];
1834  Float_t fracleri = fTCardCorrInduceEnerFrac[2][imod]+amp*fTCardCorrInduceEnerFracP1[2][imod];
1835  Float_t frac2nd = fTCardCorrInduceEnerFrac[3][imod]+amp*fTCardCorrInduceEnerFracP1[3][imod];
1836 
1837  AliDebug(1,Form("Fraction for SM %d:\n\t up-down : c %2.2e, p1 %2.2e, p2 %2.2f, sig %2.2f, fraction %2.3f\n"
1838  "\t up-down-lr: c %2.2e, p1 %2.2e, p2 %2.2f, sig %2.2f, fraction %2.3f\n"
1839  "\t left-right: c %2.2e, p1 %2.2e, p2 %2.2f, sig %2.2f, fraction %2.3f\n"
1840  "\t 2nd row : c %2.2e, p1 %2.2e, p2 %2.2f, sig %2.2f, fraction %2.3f", imod,
1841  fTCardCorrInduceEner[0][imod],fTCardCorrInduceEnerFrac[0][imod],fTCardCorrInduceEnerFracP1[0][imod],fTCardCorrInduceEnerFracWidth[0][imod],fracupdown,
1842  fTCardCorrInduceEner[1][imod],fTCardCorrInduceEnerFrac[1][imod],fTCardCorrInduceEnerFracP1[1][imod],fTCardCorrInduceEnerFracWidth[1][imod],fracupdownleri,
1843  fTCardCorrInduceEner[2][imod],fTCardCorrInduceEnerFrac[2][imod],fTCardCorrInduceEnerFracP1[2][imod],fTCardCorrInduceEnerFracWidth[2][imod],fracleri,
1844  fTCardCorrInduceEner[3][imod],fTCardCorrInduceEnerFrac[3][imod],fTCardCorrInduceEnerFracP1[3][imod],fTCardCorrInduceEnerFracWidth[3][imod],frac2nd));
1845 
1846  // Randomize the induced fraction, if requested
1847  if(fRandomizeTCard)
1848  {
1849  fracupdown = fRandom.Gaus(fracupdown ,fTCardCorrInduceEnerFracWidth[0][imod]);
1850  fracupdownleri = fRandom.Gaus(fracupdownleri,fTCardCorrInduceEnerFracWidth[1][imod]);
1851  fracleri = fRandom.Gaus(fracleri ,fTCardCorrInduceEnerFracWidth[2][imod]);
1852  frac2nd = fRandom.Gaus(frac2nd ,fTCardCorrInduceEnerFracWidth[3][imod]);
1853 
1854  AliDebug(1,Form("Randomized fraction: up-down %2.3f; up-down-left-right %2.3f; left-right %2.3f; 2nd row %2.3f",
1855  fracupdown,fracupdownleri,fracleri,frac2nd));
1856  }
1857 
1858  // Calculate induced energy
1859  Float_t indEupdown = fTCardCorrInduceEner[0][imod]+amp*fracupdown;
1860  Float_t indEupdownleri = fTCardCorrInduceEner[1][imod]+amp*fracupdownleri;
1861  Float_t indEleri = fTCardCorrInduceEner[2][imod]+amp*fracleri;
1862  Float_t indE2nd = fTCardCorrInduceEner[3][imod]+amp*frac2nd;
1863 
1864  AliDebug(1,Form("Induced energy: up-down %2.3f; up-down-left-right %2.3f; left-right %2.3f; 2nd row %2.3f",
1865  indEupdown,indEupdownleri,indEleri,indE2nd));
1866 
1867  // Check if we induce too much energy, in such case use a constant value
1868  if ( fTCardCorrMaxInduced < indE2nd ) indE2nd = fTCardCorrMaxInduced;
1869  if ( fTCardCorrMaxInduced < indEupdownleri ) indEupdownleri = fTCardCorrMaxInduced;
1870  if ( fTCardCorrMaxInduced < indEupdown ) indEupdown = fTCardCorrMaxInduced;
1871  if ( fTCardCorrMaxInduced < indEleri ) indEleri = fTCardCorrMaxInduced;
1872 
1873  AliDebug(1,Form("Induced energy, saturated?: up-down %2.3f; up-down-left-right %2.3f; left-right %2.3f; 2nd row %2.3f",
1874  indEupdown,indEupdownleri,indEleri,indE2nd));
1875 
1876  //
1877  // Add the induced energy, check if cell existed
1878  if ( okup )
1879  {
1880  fTCardCorrCellsEner[absIDup] += indEupdown;
1881 
1882  if ( fCaloCells->GetCellAmplitude(absIDup) < 0.01 ) fTCardCorrCellsNew[absIDup] = kTRUE;
1883  }
1884 
1885  if ( okdo )
1886  {
1887  fTCardCorrCellsEner[absIDdo] += indEupdown;
1888 
1889  if ( fCaloCells->GetCellAmplitude(absIDdo) < 0.01 ) fTCardCorrCellsNew[absIDdo] = kTRUE;
1890  }
1891 
1892  if ( oklr )
1893  {
1894  fTCardCorrCellsEner[absIDlr] += indEleri;
1895 
1896  if ( fCaloCells->GetCellAmplitude(absIDlr) < 0.01 ) fTCardCorrCellsNew[absIDlr] = kTRUE;
1897  }
1898 
1899  if ( okuplr )
1900  {
1901  fTCardCorrCellsEner[absIDuplr] += indEupdownleri;
1902 
1903  if ( fCaloCells->GetCellAmplitude(absIDuplr ) < 0.01 ) fTCardCorrCellsNew[absIDuplr] = kTRUE;
1904  }
1905 
1906  if ( okdolr )
1907  {
1908  fTCardCorrCellsEner[absIDdolr] += indEupdownleri;
1909 
1910  if ( fCaloCells->GetCellAmplitude(absIDdolr ) < 0.01 ) fTCardCorrCellsNew[absIDdolr] = kTRUE;
1911  }
1912 
1913  if ( okup2 )
1914  {
1915  fTCardCorrCellsEner[absIDup2] += indE2nd;
1916 
1917  if ( fCaloCells->GetCellAmplitude(absIDup2) < 0.01 ) fTCardCorrCellsNew[absIDup2] = kTRUE;
1918  }
1919 
1920  if ( okup2lr )
1921  {
1922  fTCardCorrCellsEner[absIDup2lr] += indE2nd;
1923 
1924  if ( fCaloCells->GetCellAmplitude(absIDup2lr) < 0.01 ) fTCardCorrCellsNew[absIDup2lr] = kTRUE;
1925  }
1926 
1927  if ( okdo2 )
1928  {
1929  fTCardCorrCellsEner[absIDdo2] += indE2nd;
1930 
1931  if ( fCaloCells->GetCellAmplitude(absIDdo2) < 0.01 ) fTCardCorrCellsNew[absIDdo2] = kTRUE;
1932  }
1933 
1934  if ( okdo2lr )
1935  {
1936  fTCardCorrCellsEner[absIDdo2lr] += indE2nd;
1937 
1938  if ( fCaloCells->GetCellAmplitude(absIDdo2lr) < 0.01 ) fTCardCorrCellsNew[absIDdo2lr] = kTRUE;
1939  }
1940 
1941  //
1942  // Subtract the added energy to main cell, if energy conservation is requested
1944  {
1945  if ( oklr ) fTCardCorrCellsEner[id] -= indEleri;
1946  if ( okuplr ) fTCardCorrCellsEner[id] -= indEupdownleri;
1947  if ( okdolr ) fTCardCorrCellsEner[id] -= indEupdownleri;
1948  if ( okup ) fTCardCorrCellsEner[id] -= indEupdown;
1949  if ( okdo ) fTCardCorrCellsEner[id] -= indEupdown;
1950  if ( okup2 ) fTCardCorrCellsEner[id] -= indE2nd;
1951  if ( okup2lr ) fTCardCorrCellsEner[id] -= indE2nd;
1952  if ( okdo2 ) fTCardCorrCellsEner[id] -= indE2nd;
1953  if ( okdo2lr ) fTCardCorrCellsEner[id] -= indE2nd;
1954  } // conserve energy
1955 
1956  } // cell loop
1957 
1958 }
1959 
1960 //_______________________________________________________
1962 //_______________________________________________________
1964 {
1965  AliInfo(Form("Geometry: name <%s>, matrix set <%d>, load matrix <%d>, import geo <%d> from path <%s>",
1967 
1968  if ( fAccessOCDB ) AliInfo(Form("OCDB path name <%s>", fOCDBpath.Data()));
1969  if ( fAccessOADB ) AliInfo(Form("OADB path name <%s>", fOADBFilePath.Data()));
1970 
1971  AliInfo(Form("Just Unfold clusters <%d>, new clusters list name <%s>, new cells name <%s>",
1973 
1974  if ( fFillAODFile ) AliInfo(Form("Fill new AOD file with: header <%d>, cells <%d>",fFillAODHeader,fFillAODCaloCells));
1975 
1976  AliInfo(Form("Use cell time for cluster <%d>, Apply constant time shift <%2.2f>, Do track-matching <%d>, Update cells <%d>, Input from ESD filter <%d>",
1978 
1979  AliInfo(Form("Reject events: larger than <%d>, LED <%d>, exotics <%d>", fMaxEvent, fRemoveLEDEvents, fRemoveExoticEvents));
1980 
1981  if (fCentralityBin[0] != -1 && fCentralityBin[1] != -1 )
1982  AliInfo(Form("Centrality bin [%2.2f,%2.2f], class <%s>, use AliCentrality? <%d>",
1984 
1985  if ( fSelectEMCALEvent )
1986  AliInfo(Form("Select events with signal in EMCal: E min <%2.2f>, n cell min <%d>", fEMCALEnergyCut, fEMCALNcellsCut));
1987 
1988  AliInfo(Form("MC label from cluster <%d>, Use EdepFrac <%d>, remap AODs <%d>",
1990 }
1991 
1992 //_______________________________________________________
1994 //_______________________________________________________
1996 {
1997  if(!fTCardCorrEmulation)
1998  {
1999  AliInfo("T-Card emulation not activated");
2000  return;
2001  }
2002 
2003  AliInfo(Form("T-Card emulation activated, energy conservation <%d>, randomize E <%d>, induced energy parameters:",
2005 
2006  AliInfo("T-Card emulation super-modules fraction:");
2007 
2008  for(Int_t ism = 0; ism < 22; ism++)
2009  {
2010  printf("\t sm %d, fraction %2.2f\n",ism, fTCardCorrInduceEnerProb[ism]);
2011  for(Int_t icell = 0; icell < 4; icell++)
2012  {
2013  printf("\t \t cell type %d, c %2.2e, p0 %2.2e, p1 %2.2e, sigma %2.2e \n",
2014  icell,fTCardCorrInduceEner[icell][ism],fTCardCorrInduceEnerFrac[icell][ism],fTCardCorrInduceEnerFracP1[icell][ism],fTCardCorrInduceEnerFracWidth[icell][ism]);
2015  }
2016  }
2017 }
2018 
2019 //_______________________________________________________
2023 //_______________________________________________________
2025 {
2026  Int_t j = 0;
2027  for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
2028  {
2029  AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
2030 
2031  const Int_t ncells = recPoint->GetMultiplicity();
2032  Int_t ncellsTrue = 0;
2033 
2034  if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
2035 
2036  // cells and their amplitude fractions
2037  UShort_t absIds[ncells];
2038  Double32_t ratios[ncells];
2039 
2040  //For later check embedding
2041  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2042 
2043  Float_t clusterE = 0;
2044  for (Int_t c = 0; c < ncells; c++)
2045  {
2046  AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
2047 
2048  absIds[ncellsTrue] = digit->GetId();
2049  ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
2050 
2051  if ( !fRecParam->GetUnfold() && (ratios[ncellsTrue] > 1 || ratios[ncellsTrue] < 1) )
2052  AliWarning(Form("recpoint cell E %2.3f but digit E %2.3f and no unfolding", recPoint->GetEnergiesList()[c], digit->GetAmplitude()));
2053 
2054  // In case of unfolding, remove digits with unfolded energy too low
2055  if(fSelectCell)
2056  {
2057  if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
2058  {
2059 
2060  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
2061  recPoint->GetEnergiesList()[c],digit->GetAmplitude()));
2062  continue;
2063  } // if cuts
2064  }// Select cells
2065 
2066  //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
2067  clusterE +=recPoint->GetEnergiesList()[c];
2068 
2069  // In case of embedding, fill ratio with amount of signal,
2070  if (aodIH && aodIH->GetMergeEvents())
2071  {
2072  //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
2073  AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
2074  AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
2075 
2076  Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2077  //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2078  Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2079  //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
2080 
2081  if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
2082  //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
2083 
2084  }//Embedding
2085 
2086  ncellsTrue++;
2087 
2088  }// cluster cell loop
2089 
2090  if (ncellsTrue < 1)
2091  {
2092  AliDebug(2,Form("Skipping cluster with no cells avobe threshold E = %f, ncells %d",
2093  recPoint->GetEnergy(), ncells));
2094  continue;
2095  }
2096 
2097  //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
2098 
2099  if(clusterE < fRecParam->GetClusteringThreshold())
2100  {
2101  AliDebug(2,Form("Remove cluster with energy below seed threshold %f",clusterE));
2102  continue;
2103  }
2104 
2105  TVector3 gpos;
2106  Float_t g[3];
2107 
2108  // calculate new cluster position
2109 
2110  recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
2111  recPoint->GetGlobalPosition(gpos);
2112  gpos.GetXYZ(g);
2113 
2114  // create a new cluster
2115 
2116  (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
2117  AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
2118  j++;
2119  clus->SetType(AliVCluster::kEMCALClusterv1);
2120  clus->SetE(clusterE);
2121  clus->SetPosition(g);
2122  clus->SetNCells(ncellsTrue);
2123  clus->SetCellsAbsId(absIds);
2124  clus->SetCellsAmplitudeFraction(ratios);
2125  clus->SetChi2(-1); //not yet implemented
2126  clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
2127  clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
2128  clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
2129 
2130  if(ncells == ncellsTrue)
2131  {
2132  Float_t elipAxis[2];
2133  recPoint->GetElipsAxis(elipAxis);
2134  clus->SetM02(elipAxis[0]*elipAxis[0]) ;
2135  clus->SetM20(elipAxis[1]*elipAxis[1]) ;
2136  clus->SetDispersion(recPoint->GetDispersion());
2137  }
2138  else if(fSelectCell)
2139  {
2140  // In case some cells rejected, in unfolding case, recalculate
2141  // shower shape parameters and position
2142  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
2143 
2144  AliVCaloCells* cells = 0x0;
2145  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
2146  else cells = InputEvent()->GetEMCALCells();
2147 
2148  fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
2149  fRecoUtils->RecalculateClusterPID(clus);
2150  fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
2151  }
2152 
2153  // MC
2154 
2157  else
2158  {
2159  //
2160  // Normal case, trust what the clusterizer has found
2161  //
2162  Int_t parentMult = 0;
2163  Int_t *parentList = recPoint->GetParents(parentMult);
2164  Float_t *parentListDE = recPoint->GetParentsDE(); // deposited energy
2165 
2166  clus->SetLabel(parentList, parentMult);
2168  clus->SetClusterMCEdepFractionFromEdepArray(parentListDE);
2169 
2170  //
2171  // Set the cell energy deposition fraction map:
2172  //
2173  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
2174  {
2175  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
2176 
2177  // Get the digit that originated this cell cluster
2178 // AliVCaloCells* cells = 0x0;
2179 // if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
2180 // else cells = InputEvent()->GetEMCALCells();
2181 
2182  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
2183  {
2184  Int_t idigit = fCellLabels[absIds[icell]];
2185 
2186  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
2187 
2188  // Find the 4 MC labels that contributed to the cluster and their
2189  // deposited energy in the current digit
2190 
2191  mcEdepFracPerCell[icell] = 0; // init
2192 
2193  Int_t nparents = dig->GetNiparent();
2194  if ( nparents > 0 )
2195  {
2196  Int_t digLabel =-1 ;
2197  Float_t edep = 0 ;
2198  Float_t edepTot = 0 ;
2199  Float_t mcEDepFrac[4] = {0,0,0,0};
2200 
2201  // all parents in digit
2202  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2203  {
2204  digLabel = dig->GetIparent (jndex+1);
2205  edep = dig->GetDEParent(jndex+1);
2206  edepTot += edep;
2207 
2208  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
2209  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
2210  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
2211  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
2212  } // all prarents in digit
2213 
2214  // Divide energy deposit by total deposited energy
2215  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
2216  if(edepTot > 0.01)
2217  {
2218  mcEdepFracPerCell[icell] = clus->PackMCEdepFraction(mcEDepFrac);
2219  }
2220  } // at least one parent label in digit
2221  } // cell in cluster loop
2222 
2223  clus->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
2224 
2225  delete [] mcEdepFracPerCell;
2226 
2227  } // at least one parent in cluster, do the cell primary packing
2228  }
2229  } // recPoints loop
2230 }
2231 
2232 //_____________________________________________________________________
2235 //_____________________________________________________________________
2237 {
2238  if(label < 0) return ;
2239 
2240  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
2241  if(!evt) return ;
2242 
2243  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
2244  if(!arr) return ;
2245 
2246  if(label < arr->GetEntriesFast())
2247  {
2248  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
2249  if(!particle) return ;
2250 
2251  if(label == particle->Label()) return ; // label already OK
2252  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
2253  }
2254  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
2255 
2256  // loop on the particles list and check if there is one with the same label
2257  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
2258  {
2259  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
2260  if(!particle) continue ;
2261 
2262  if(label == particle->Label())
2263  {
2264  label = ind;
2265  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
2266  return;
2267  }
2268  }
2269 
2270  label = -1;
2271 
2272  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
2273 }
2274 
2275 //________________________________________________
2277 //________________________________________________
2279 {
2280  for(Int_t j = 0; j < fgkNEMCalCells; j++)
2281  {
2282  fOrgClusterCellId[j] = -1;
2283  fCellLabels[j] = -1;
2284  fCellSecondLabels[j] = -1;
2285  fCellTime[j] = 0.;
2286  fCellMatchdEta[j] = -999;
2287  fCellMatchdPhi[j] = -999;
2288 
2289  fTCardCorrCellsEner[j] = 0.;
2290  fTCardCorrCellsNew [j] = kFALSE;
2291  }
2292 }
2293 
2294 //_____________________________________________________________________________________________________
2298 //_____________________________________________________________________________________________________
2300  AliAODCaloCluster * clus)
2301 {
2302  Int_t parentMult = 0;
2303  Int_t *parentList = recPoint->GetParents(parentMult);
2304  clus->SetLabel(parentList, parentMult);
2305 
2306  //Write the second major contributor to each MC cluster.
2307  Int_t iNewLabel ;
2308  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2309  {
2310 
2311  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2312  if(idCell>=0)
2313  {
2314  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
2315  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
2316  {
2317  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
2318  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
2319  }
2320  if (iNewLabel == 1)
2321  {
2322  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
2323  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
2324  {
2325  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
2326  }
2327 
2328  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
2329  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
2330  delete [] newLabelArray;
2331  }
2332  } // positive cell number
2333  }
2334 }
2335 
2336 //___________________________________________________________________________________________________
2340 //___________________________________________________________________________________________________
2342 {
2343  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
2344  clArray.Reset();
2345  Int_t nClu = 0;
2346  Int_t nLabTotOrg = 0;
2347  Float_t emax = -1;
2348  Int_t idMax = -1;
2349 
2350  AliVEvent * event = fEvent;
2351 
2352  // In case of embedding the MC signal is in the added event
2353  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2354  if(aodIH && aodIH->GetEventToMerge()) //Embedding
2355  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
2356 
2357 
2358  // Find the clusters that originally had the cells
2359  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2360  {
2361  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2362 
2363  if(idCell>=0)
2364  {
2365  Int_t idCluster = fOrgClusterCellId[idCell];
2366 
2367  Bool_t set = kTRUE;
2368  for(Int_t icl =0; icl < nClu; icl++)
2369  {
2370  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
2371  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
2372  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
2373  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
2374  }
2375  if( set && idCluster >= 0)
2376  {
2377  clArray.SetAt(idCluster,nClu++);
2378  //printf("******** idCluster %d \n",idCluster);
2379  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
2380 
2381  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
2382 
2383  //Search highest E cluster
2384  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2385  //printf("\t E %f\n",clOrg->E());
2386  if(emax < clOrg->E())
2387  {
2388  emax = clOrg->E();
2389  idMax = idCluster;
2390  }
2391  }
2392  }
2393  }// cell loop
2394 
2395  if(nClu==0 || nLabTotOrg == 0)
2396  {
2397  //if(clus->E() > 0.25) printf("AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters() - Check: N org clusters %d, n tot labels %d, cluster E %f, n cells %d\n",nClu,nLabTotOrg,clus->E(), clus->GetNCells());
2398  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
2399  //printf("\n");
2400  }
2401 
2402  // Put the first in the list the cluster with highest energy
2403  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
2404  {
2405  Int_t maxIndex = -1;
2406  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
2407  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2408  {
2409  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
2410  }
2411 
2412  if(firstCluster >=0 && idMax >=0)
2413  {
2414  clArray.SetAt(idMax,0);
2415  clArray.SetAt(firstCluster,maxIndex);
2416  }
2417  }
2418 
2419  // Get the labels list in the original clusters, assign all to the new cluster
2420  TArrayI clMCArray(nLabTotOrg) ;
2421  clMCArray.Reset();
2422 
2423  Int_t nLabTot = 0;
2424  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2425  {
2426  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
2427  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
2428  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2429  Int_t nLab = clOrg->GetNLabels();
2430 
2431  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
2432  {
2433  Int_t lab = clOrg->GetLabelAt(iLab) ;
2434  if(lab>=0)
2435  {
2436  Bool_t set = kTRUE;
2437  //printf("\t \t Set Label %d \n", lab);
2438  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
2439  {
2440  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
2441  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
2442  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
2443  }
2444  if( set ) clMCArray.SetAt(lab,nLabTot++);
2445  }
2446  }
2447  }// cluster loop
2448 
2449  // Set the final list of labels
2450 
2451  clus->SetLabel(clMCArray.GetArray(), nLabTot);
2452 
2453 // printf("Final list of labels for new cluster : \n");
2454 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
2455 // {
2456 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
2457 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
2458 // printf(" org %d ",label);
2459 // RemapMCLabelForAODs(label);
2460 // printf(" new %d \n",label);
2461 // }
2462 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
2463 }
2464 
2465 //____________________________________________________________
2466 // Init geometry, create list of output clusters and cells.
2467 //____________________________________________________________
2469 {
2470  // Clusters
2471  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2472 
2473  if(fOutputAODBranchName.Length()==0)
2474  {
2475  fOutputAODBranchName = "newEMCALClustersArray";
2476  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
2477  }
2478 
2480 
2481  if( fFillAODFile )
2482  {
2483  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2484 
2485  //fOutputAODBranch->SetOwner(kFALSE);
2486 
2487  AddAODBranch("TClonesArray", &fOutputAODBranch);
2488  }
2489 
2490  // Cells
2491  if ( fOutputAODBranchName.Length() > 0 )
2492  {
2493  fOutputAODCells = new AliAODCaloCells(fOutputAODCellsName,fOutputAODCellsName,AliAODCaloCells::kEMCALCell);
2494 
2495  if( fFillAODFile ) AddAODBranch("AliAODCaloCells", &fOutputAODCells);
2496  }
2497 }
2498 
2499 //_______________________________________________________________________
2502 //________________________________________________________________________
2504 {
2505  // Update cells only in case re-calibration was done
2506  // or bad map applied or additional T-Card cells added.
2507  if(!fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
2508  !fRecoUtils->IsRecalibrationOn() &&
2509  !fRecoUtils->IsRunDepRecalibrationOn() &&
2510  !fRecoUtils->IsTimeRecalibrationOn() &&
2511  !fRecoUtils->IsL1PhaseInTimeRecalibrationOn() &&
2512  !fTCardCorrEmulation ) return;
2513 
2514  const Int_t ncells = fCaloCells->GetNumberOfCells();
2515  const Int_t ndigis = fDigitsArr->GetEntries();
2516 
2517  if ( fOutputAODCellsName.Length() > 0 )
2518  {
2519  fOutputAODCells->DeleteContainer();
2520  fOutputAODCells->CreateContainer(ndigis);
2521  }
2522  else if ( ncells != ndigis ) // update case
2523  {
2524  fCaloCells->DeleteContainer();
2525  fCaloCells->CreateContainer(ndigis);
2526  }
2527 
2528  for (Int_t idigit = 0; idigit < ndigis; ++idigit)
2529  {
2530  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
2531 
2532  Double_t cellAmplitude = digit->GetAmplitude();
2533  Short_t cellNumber = digit->GetId();
2534  Double_t cellTime = digit->GetTime();
2535 
2536  Bool_t highGain = kFALSE;
2537  if( digit->GetType() == AliEMCALDigit::kHG ) highGain = kTRUE;
2538 
2539  // Only for MC
2540  // Get the label of the primary particle that generated the cell
2541  // Assign the particle that deposited more energy
2542  Int_t nparents = digit->GetNiparent();
2543  Int_t cellMcEDepFrac =-1 ;
2544  Float_t cellMcLabel =-1.;
2545  if ( nparents > 0 )
2546  {
2547  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2548  {
2549  if(cellMcEDepFrac >= digit->GetDEParent(jndex+1)) continue ;
2550 
2551  cellMcLabel = digit->GetIparent (jndex+1);
2552  cellMcEDepFrac= digit->GetDEParent(jndex+1);
2553  } // all primaries in digit
2554  } // select primary label
2555 
2556  if ( cellMcEDepFrac < 0 ) cellMcEDepFrac = 0.;
2557 
2558  if ( fUpdateCell )
2559  fCaloCells ->SetCell(idigit, cellNumber, cellAmplitude, cellTime, cellMcLabel, cellMcEDepFrac, highGain);
2560  else
2561  fOutputAODCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime, cellMcLabel, cellMcEDepFrac, highGain);
2562  }
2563 }
2564 
2565 //_______________________________________________________
2576 //_______________________________________________________
2578 {
2579  //-------
2580  // Step 1
2581 
2582  // Remove the contents of AOD branch output list set in the previous event
2583  fOutputAODBranch->Clear("C");
2584 
2585  LoadBranches();
2586 
2587  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
2588  // If we need a centrality bin, we select only those events in the corresponding bin.
2589  if( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
2590  {
2591  Float_t cen = GetEventCentrality();
2592  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
2593  }
2594 
2595  // Intermediate array with new clusters : init the array only once or clear from previous event
2596  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
2597  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
2598 
2599 
2600  // In case of analysis in multiple runs, check the OADB again
2601  if ( InputEvent()->GetRunNumber() != fRun )
2602  {
2603  fRun = InputEvent()->GetRunNumber();
2604 
2605  fOADBSet = kFALSE; // recover the OADB for this run
2606 
2607  AliInfo(Form("Set run to %d",fRun));
2608  }
2609 
2610  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
2611 
2612  // Get the event, do some checks and settings
2613  CheckAndGetEvent() ;
2614 
2615  if (!fEvent)
2616  {
2617  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
2618  return ;
2619  }
2620 
2621  // Init pointers, geometry, clusterizer, ocdb, aodb
2622 
2623  if(fAccessOCDB) AccessOCDB();
2624  if(fAccessOADB) AccessOADB(); // only once
2625 
2627 
2628  // Print once the analysis parameters
2629  if ( fDebug > 0 || !fPrintOnce )
2630  {
2631  fRecParam->Print("reco");
2632 
2633  PrintParam();
2634 
2635  PrintTCardParam();
2636 
2637  fPrintOnce = kTRUE;
2638  }
2639 
2640  //-------
2641  // Step 2
2642 
2643  // Make clusters
2645  else ClusterizeCells() ;
2646 
2647  //-------
2648  // Step 3
2649 
2651 
2652  if ( fUpdateCell || fOutputAODCellsName.Length() > 0 )
2653  UpdateCells();
2654 
2655 }
2656 
2657 
TString fOutputAODCellsName
New of output cells AOD branch.
void FillAODCaloCells()
Put calo cells in standard branch.
Float_t fTCardCorrMinAmp
Minimum cell energy to induce signal on adjacent cells.
AliVCaloCells * fCaloCells
! CaloCells container
Bool_t fTCardCorrEmulation
Activate T-Card cells energy correlation.
Float_t fTCardCorrInduceEnerFracP1[4][22]
Induced energy loss gauss fraction param1 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
double Double_t
Definition: External.C:58
virtual void ResetArrays()
Reset arrays containing information for all possible cells.
Definition: External.C:236
const int debug
Definition: scanAll.C:15
TString GetPass()
Get or guess pass number/string from path of filename.
Bool_t fFillAODHeader
Copy header to standard branch.
TClonesArray * fDigitsArr
! Digits array
Bool_t fTCardCorrClusEnerConserv
When making correlation, subtract from the reference cell the induced energy on the neighbour cells...
AliEMCALClusterizer * fClusterizer
! EMCAL clusterizer
Bool_t fJustUnfold
Just unfold, do not recluster.
Bool_t fImportGeometryFromFile
Import geometry settings in geometry.root file.
Bool_t fSelectCell
Reject cells from cluster if energy is too low and recalculate position/energy and other...
TString fConfigName
Name of analysis configuration file.
TCanvas * c
Definition: TestFitELoss.C:172
void PrintTCardParam()
Print parameters for T-Card correlation emulation.
Int_t fCellSecondLabels[fgkNEMCalCells]
Array with Second MC label to be passed to digit.
void SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint *recPoint, AliAODCaloCluster *clus)
Float_t fCellMatchdEta[fgkNEMCalCells]
Array with cluster-track dPhi.
TClonesArray * fOutputAODBranch
! AOD Branch with output clusters
TString fOADBFilePath
Default path $ALICE_PHYSICS/OADB/EMCAL, if needed change.
Float_t fTCardCorrCellsEner[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells.
TString fImportGeometryFilePath
path fo geometry.root file
Float_t fTCardCorrInduceEnerFrac[4][22]
Induced energy loss gauss fraction param0 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
void ConfigureEMCALRecoUtils(Bool_t bMC=kFALSE, Bool_t bExotic=kTRUE, Bool_t bNonLin=kFALSE, Bool_t bRecalE=kTRUE, Bool_t bBad=kTRUE, Bool_t bRecalT=kTRUE, Int_t debug=-1)
Bool_t fPrintOnce
Print once analysis parameters.
Bool_t fSelectEMCALEvent
Process the event if there is some high energy cluster.
Float_t fTCardCorrMaxInduced
Maximum induced energy signal on adjacent cells.
Float_t GetEventCentrality() const
Get centrality/multiplicity percentile.
Bool_t fInputFromFilter
Get the input from AODs from the filter.
TString fOutputAODBranchName
New of output clusters AOD branch.
AliEMCALGeometry * fGeom
EMCAL geometry.
void ClusterUnfolding()
Take the event clusters and unfold them.
int Int_t
Definition: External.C:63
Bool_t fAccessOADB
Get calibration from OADB for EMCAL.
Bool_t fRemoveLEDEvents
Remove LED events, use only for LHC11a.
unsigned int UInt_t
Definition: External.C:33
AliEMCALRecoUtils * fRecoUtils
Access to factorized reconstruction algorithms.
float Float_t
Definition: External.C:68
Bool_t fRejectBelowThreshold
split (false-default) or reject (true) cell energy below threshold after UF
Float_t fConstantTimeShift
Apply a 600 ns time shift in case of simulation, shift in ns.
void SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster *clus)
Bool_t fUpdateCell
On/Off the upate of the CaloCells container.
Float_t fSelectCellMinE
Min energy cell threshold, after unfolding.
Bool_t fOADBSet
AODB parameters already set.
AliAODCaloCells * fOutputAODCells
! AOD Branch with output cells
Float_t fTCardCorrInduceEner[4][22]
Induced energy loss gauss constant on 0-same row, diff col, 1-up/down cells left/right col 2-left/rig...
Int_t fCellLabels[fgkNEMCalCells]
Array with MC label to be passed to digit.
Bool_t fRemoveExoticEvents
Remove exotic events.
Bool_t fDoTrackMatching
On/Off the matching recalculation to speed up analysis in PbPb.
Float_t fCellMatchdPhi[fgkNEMCalCells]
Array with cluster-track dEta.
Float_t fTCardCorrInduceEnerProb[22]
Probability to induce energy loss per SM.
TGeoHMatrix * fGeomMatrix[22]
Geometry matrices with alignments.
Bool_t IsLEDEvent(const Int_t run)
Check if event is LED, is so remove it. Affected LHC11a runs.
AliEMCALRecParam * fRecParam
Reconstruction parameters container.
TString fGeomName
Name of geometry to use.
Bool_t fOutputAODBranchSet
Set the AOD clusters branch in the input event once.
void PrintParam()
Print clusterization task parameters.
short Short_t
Definition: External.C:23
TString fOCDBpath
Path with OCDB location.
AliEMCALAfterBurnerUF * fUnfolder
! Unfolding procedure
Bool_t fGeomMatrixSet
Set geometry matrices only once, for the first event.
Float_t fCentralityBin[2]
Minimum and maximum value of the centrality for the analysis.
Bool_t fFillAODCaloCells
Copy calocells to standard branch.
TObjArray * fCaloClusterArr
! CaloClusters array
Double_t fCellTime[fgkNEMCalCells]
Array with cluster time to be passed to digit in case of AODs.
TString fCentralityClass
Name of selected centrality class.
TFile * file
TList with histograms for a given trigger.
Bool_t fAccessOCDB
Need to access info from OCDB (not really)
TObjArray * fClusterArr
! Recpoints array
unsigned short UShort_t
Definition: External.C:28
Bool_t fRecalibrateWithClusterTime
Use fCellTime to store time of cells in cluster.
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void FillAODHeader()
Put event header information in standard AOD branch.
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)...
Int_t fEMCALNcellsCut
At least an EMCAL cluster with fNCellsCut cells over fEnergyCut.
Reclusterize EMCal clusters, put them in a new branch for other following analysis.
Bool_t fTCardCorrCellsNew[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells, that before had no signal.
Bool_t fLoadGeomMatrices
Matrices set from configuration, not get from geometry.root or from ESDs/AODs.
Float_t fEMCALEnergyCut
At least an EMCAL cluster with this energy in the event.
Bool_t fRandomizeTCard
Use random induced energy.
Bool_t fUseAliCentrality
Use the centrality estimator from AliCentrality or AliMultSelection.
Float_t fTCardCorrInduceEnerFracWidth[4][22]
Induced energy loss gauss witdth on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ ...
Int_t fOrgClusterCellId[fgkNEMCalCells]
Array ID of cluster to wich the cell belongs in unmodified clusters.
Float_t fSelectCellMinFrac
Min fraction of cell energy after unfolding cut.
Bool_t fRemapMCLabelForAODs
Remap AOD cells MC label. Needed in old AOD productions.