AliPhysics  5364b50 (5364b50)
AliCalorimeterUtils.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 system ---
17 #include <TGeoManager.h>
18 #include <TH2F.h>
19 #include <TCanvas.h>
20 #include <TStyle.h>
21 #include <TPad.h>
22 #include <TFile.h>
23 #include <TParticle.h>
24 #include <AliMCEvent.h>
25 
26 // --- ANALYSIS system ---
27 #include "AliCalorimeterUtils.h"
28 #include "AliDataFile.h"
29 #include "AliESDEvent.h"
30 #include "AliMCEvent.h"
31 #include "AliCaloTrackParticle.h"
32 #include "AliVCluster.h"
33 #include "AliVCaloCells.h"
34 #include "AliAODCaloCluster.h"
35 #include "AliOADBContainer.h"
36 #include "AliAnalysisManager.h"
37 #include "AliAODMCParticle.h"
38 #include "AliVParticle.h"
39 #include "AliLog.h"
40 
41 // --- Detector ---
42 #include "AliEMCALGeometry.h"
43 #include "AliPHOSGeoUtils.h"
44 
45 #include "AliFiducialCut.h" // Needed for detector flag enum kEMCAL, kPHOS
46 
48 ClassImp(AliCalorimeterUtils) ;
50 
51 
52 //____________________________________________
54 //____________________________________________
56 TObject(), fDebug(0),
57 fEMCALGeoName(""),
58 fPHOSGeoName (""),
59 fEMCALGeo(0x0), fPHOSGeo(0x0),
60 fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
61 fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
62 fRemoveBadChannels(kFALSE), fPHOSBadChannelMap(0x0),
63 fNCellsFromPHOSBorder(0),
64 fNMaskCellColumns(0), fMaskCellColumns(0x0),
65 fRecalibration(kFALSE), fRunDependentCorrection(kFALSE),
66 fPHOSRecalibrationFactors(), fEMCALRecoUtils(new AliEMCALRecoUtils),
67 fRecalculatePosition(kFALSE), fCorrectELinearity(kFALSE),
68 fRecalculateMatching(kFALSE),
69 fCutR(20), fCutZ(20),
70 fCutEta(20), fCutPhi(20),
71 fLocMaxCutE(0), fLocMaxCutEDiff(0),
72 fPlotCluster(0), fOADBSet(kFALSE),
73 fOADBForEMCAL(kFALSE), fOADBForPHOS(kFALSE),
74 fOADBFilePathEMCAL(""), fOADBFilePathPHOS(""),
75 fImportGeometryFromFile(0), fImportGeometryFilePath(""),
76 fNSuperModulesUsed(0),
77 fFirstSuperModuleUsed(-1), fLastSuperModuleUsed(-1),
78 fRunNumber(0),
79 fMCECellClusFracCorrOn(0), fMCECellClusFracCorrParam()
80 {
82  for(Int_t i = 0; i < 22; i++) fEMCALMatrix[i] = 0 ;
83  for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
84 }
85 
86 //_________________________________________
87 // Destructor.
88 //_________________________________________
90 {
91  //if(fPHOSGeo) delete fPHOSGeo ;
92  if(fEMCALGeo) delete fEMCALGeo ;
93 
95 {
96  fPHOSBadChannelMap->Clear();
97  delete fPHOSBadChannelMap;
98  }
99 
101 {
102  fPHOSRecalibrationFactors->Clear();
104  }
105 
106  if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
107  if(fNMaskCellColumns) delete [] fMaskCellColumns;
108 }
109 
110 //____________________________________________________
113 //____________________________________________________
114 void AliCalorimeterUtils::AccessOADB(AliVEvent* event)
115 {
116  // Set it only once
117  if(fOADBSet) return ;
118 
119  if(fRunNumber <= 0) fRunNumber = event->GetRunNumber() ; // take the run number from the event itself
120  TString pass = GetPass();
121 
122  // EMCAL
123  if(fOADBForEMCAL)
124  {
125  AliInfo(Form("Get AODB parameters from EMCAL in %s for run %d, and <%s>",fOADBFilePathEMCAL.Data(),fRunNumber,pass.Data()));
126 
127  Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
128 
129  // Bad map
131  {
132  AliOADBContainer *contBC=new AliOADBContainer("");
133  if(fOADBFilePathEMCAL!="")
134  contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePathEMCAL.Data()),"AliEMCALBadChannels");
135  else
136  contBC->InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALBadChannels.root").data(),"AliEMCALBadChannels");
137 
138  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRunNumber);
139 
140  if(arrayBC)
141  {
143  AliInfo("Remove EMCAL bad cells");
144 
145  for (Int_t i=0; i<nSM; ++i)
146  {
147  TH2I *hbm = GetEMCALChannelStatusMap(i);
148 
149  if (hbm)
150  delete hbm;
151 
152  hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
153 
154  if (!hbm)
155  {
156  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
157  continue;
158  }
159 
160  hbm->SetDirectory(0);
162 
163  } // loop
164  } else AliInfo("Do NOT remove EMCAL bad channels\n"); // run array
165 
166  delete contBC;
167  } // Remove bad
168 
169  // Energy Recalibration
170  if(fRecalibration)
171  {
172  AliOADBContainer *contRF=new AliOADBContainer("");
173 
174  if(fOADBFilePathEMCAL!="")
175  contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRecalib");
176  else
177  contRF->InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALRecalib.root").data(),"AliEMCALRecalib");
178 
179  TObjArray *recal=(TObjArray*)contRF->GetObject(fRunNumber);
180 
181  if(recal)
182  {
183  TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
184 
185  if(recalpass)
186  {
187  TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
188 
189  if(recalib)
190  {
191  AliInfo("Recalibrate EMCAL");
192  for (Int_t i=0; i < nSM; ++i)
193  {
195 
196  if (h)
197  delete h;
198 
199  h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
200 
201  if (!h)
202  {
203  AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
204  continue;
205  }
206 
207  h->SetDirectory(0);
208 
210  } // SM loop
211  } else AliInfo("Do NOT recalibrate EMCAL, no params object array"); // array ok
212  } else AliInfo("Do NOT recalibrate EMCAL, no params for pass"); // array pass ok
213  } else AliInfo("Do NOT recalibrate EMCAL, no params for run"); // run number array ok
214 
215  delete contRF;
216  // once set, apply run dependent corrections if requested
217  //fEMCALRecoUtils->SetRunDependentCorrections(fRunNumber);
218 
219  } // Recalibration on
220 
221  // Energy Recalibration, apply on top of previous calibration factors
223  {
224  AliOADBContainer *contRFTD=new AliOADBContainer("");
225 
226  if(fOADBFilePathEMCAL!="")
227  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRunDepTempCalibCorrections");
228  else
229  contRFTD->InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALTemperatureCorrCalib.root").data(),"AliEMCALRunDepTempCalibCorrections");
230 
231  TH1S *htd=(TH1S*)contRFTD->GetObject(fRunNumber);
232 
233  //If it did not exist for this run, get closes one
234  if (!htd)
235  {
236  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",fRunNumber));
237  // let's get the closest fRunNumber instead then..
238  Int_t lower = 0;
239  Int_t ic = 0;
240  Int_t maxEntry = contRFTD->GetNumberOfEntries();
241 
242  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < fRunNumber) )
243  {
244  lower = ic;
245  ic++;
246  }
247 
248  Int_t closest = lower;
249  if ( (ic<maxEntry) &&
250  (contRFTD->LowerLimit(ic)-fRunNumber) < (fRunNumber - contRFTD->UpperLimit(lower)) )
251  {
252  closest = ic;
253  }
254 
255  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
256  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
257  }
258 
259  if(htd)
260  {
261  AliInfo("Recalibrate (Temperature) EMCAL");
262 
263  for (Int_t ism=0; ism<nSM; ++ism)
264  {
265  for (Int_t icol=0; icol<48; ++icol)
266  {
267  for (Int_t irow=0; irow<24; ++irow)
268  {
269  Float_t factor = GetEMCALChannelRecalibrationFactor(ism,icol,irow);
270 
271  Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
272  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
273 
274  //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
275  // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
276 
277  SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
278  } // columns
279  } // rows
280  } // SM loop
281  } else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
282 
283  delete contRFTD;
284  } // Run by Run T calibration
285 
286  // Time Recalibration
287  if(fEMCALRecoUtils->IsTimeRecalibrationOn())
288  {
289  AliOADBContainer *contTRF=new AliOADBContainer("");
290 
291  if(fOADBFilePathEMCAL!="")
292  contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeCalib");
293  else
294  contTRF->InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALTimeCalib.root").data(),"AliEMCALTimeCalib");
295 
296  TObjArray *trecal=(TObjArray*)contTRF->GetObject(fRunNumber);
297 
298  if(trecal)
299  {
300  // pass number should be pass1 except on Run1 and special cases
301  TString passM = pass;
302  if ( pass=="spc_calo" ) passM = "pass3";
303  if ( fRunNumber > 209121 ) passM = "pass1";//run2 periods
304  if ( pass == "muon_calo_pass1" && fRunNumber > 209121 && fRunNumber < 244284 )
305  passM = "pass0";//period LHC15a-m
306 
307  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
308 
309  if(trecalpass)
310  {
311  AliInfo("Time Recalibrate EMCAL");
312  for (Int_t ibc = 0; ibc < 4; ++ibc)
313  {
315 
316  if (h)
317  delete h;
318 
319  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
320 
321  if (!h)
322  {
323  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
324  continue;
325  }
326 
327  h->SetDirectory(0);
328 
330  } // bunch crossing loop
331  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
332  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
333 
334  delete contTRF;
335  } // Time Recalibration on
336 
337  // Time L1 phase racalibration
338  if(fEMCALRecoUtils->IsL1PhaseInTimeRecalibrationOn())
339  {
340  AliOADBContainer *contTRF=new AliOADBContainer("");
341  if(fOADBFilePathEMCAL!="")
342  contTRF->InitFromFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeL1PhaseCalib");
343  else
344  contTRF->InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALTimeL1PhaseCalib.root").data(),"AliEMCALTimeL1PhaseCalib");
345 
346  TObjArray *trecal=(TObjArray*)contTRF->GetObject(fRunNumber);
347  if(!trecal)
348  {
349  AliError(Form("L1 phase time recal: No params for run %d. Default used.",fRunNumber)); // run number array ok
350  trecal=(TObjArray*)contTRF->GetObject(0); // Try default object
351  }
352 
353  if(trecal)
354  {
355  // Only 1 L1 phase correction possible, except special cases
356  TString passM = "pass1";
357 
358  if ( pass == "muon_calo_pass1" && fRunNumber > 209121 && fRunNumber < 244284 )
359  passM = "pass0";//period LHC15a-m
360 
361  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
362  if(!trecalpass)
363  {
364  AliInfo(Form("L1 phase time recal: No params for run %d and pass %s, try default", fRunNumber, passM.Data()));
365 
366  trecal->Delete();
367 
368  trecal=(TObjArray*)contTRF->GetObject(0);
369 
370  if(trecal)
371  trecalpass=(TObjArray*)trecal->FindObject("pass1");
372 
373  AliInfo("Time L1 phase Recalibrate EMCAL");
374  }
375 
376  if(trecalpass)
377  {
379 
380  if (h) delete h;
381 
382  h = (TH1C*)trecalpass->FindObject(Form("h%d",fRunNumber));
383 
384  if (!h) AliError(Form("Could not load h%d",fRunNumber));
385 
386  h->SetDirectory(0);
387 
389  }
390  else
391  {
392  AliError("Do not calibrate L1 phase time");
393  fEMCALRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
394  }
395  }
396  else
397  {
398  AliError("Do not calibrate L1 phase time");
399  fEMCALRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
400  }
401 
402  delete contTRF;
403  }//End of Time L1 phase racalibration
404 
405  }// EMCAL
406 
407  // PHOS
408  if(fOADBForPHOS)
409  {
410  AliInfo(Form("Get AODB parameters from PHOS in %s for run %d, and <%s>",fOADBFilePathPHOS.Data(),fRunNumber,pass.Data()));
411 
412  // Bad map
414  {
415  AliOADBContainer badmapContainer(Form("phosBadMap"));
416  TString fileName="$ALICE_PHYSICS/OADB/PHOS/PHOSBadMaps.root";
417  badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
418 
419  //Use a fixed run number from year 2010, this year not available yet.
420  TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
421  if(!maps)
422  {
423  AliInfo(Form("Can not read PHOS bad map for run %d",fRunNumber)) ;
424  }
425  else
426  {
427  AliInfo(Form("Setting PHOS bad map with name %s",maps->GetName())) ;
428 
429  for(Int_t mod = 1; mod < 5; mod++)
430  {
431  TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
432 
433  if(hbmPH)
434  delete hbmPH ;
435 
436  hbmPH = (TH2I*)maps->At(mod);
437 
438  if(hbmPH) hbmPH->SetDirectory(0);
439 
440  SetPHOSChannelStatusMap(mod-1,hbmPH);
441 
442  } // modules loop
443  } // maps exist
444  } // Remove bad channels
445  } // PHOS
446 
447  // Parameters already set once, so do not it again
448  fOADBSet = kTRUE;
449 }
450 
451 //_____________________________________________________________
454 //_____________________________________________________________
455 void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent)
456 {
457  // First init the geometry, a priory not done before
458  if(fRunNumber <=0 ) fRunNumber = inputEvent->GetRunNumber() ;
459 
460  InitPHOSGeometry ();
462 
463  //Get the EMCAL transformation geometry matrices from ESD
465  {
467  {
468  AliInfo("Load user defined EMCAL geometry matrices");
469 
470  // OADB if available
471  AliOADBContainer emcGeoMat("AliEMCALgeo");
472  if(fOADBFilePathEMCAL!="")
473  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
474  else
475  emcGeoMat.InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALlocal2master.root").data(),"AliEMCALgeo");
476 
477  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(fRunNumber,"EmcalMatrices");
478 
479  for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
480  {
481  if (!fEMCALMatrix[mod]) // Get it from OADB
482  {
483  AliDebug(1,Form("EMCAL matrices SM %d, %p", mod,((TGeoHMatrix*) matEMCAL->At(mod))));
484  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
485 
486  fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
487  }
488 
489  if(fEMCALMatrix[mod])
490  {
491  if(fDebug > 1)
492  fEMCALMatrix[mod]->Print();
493 
494  fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
495  }
496  else if(gGeoManager)
497  {
498  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
499  fEMCALGeo->SetMisalMatrix(fEMCALGeo->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
500  }
501  else
502  {
503  AliError(Form("Alignment atrix for SM %d is not available",mod));
504  }
505  } // SM loop
506 
507  fEMCALGeoMatrixSet = kTRUE;//At least one, so good
508 
509  } // Load matrices
510  else if (!gGeoManager)
511  {
512  AliDebug(1,"Load EMCAL misalignment matrices");
513  if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
514  {
515  for(Int_t mod = 0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
516  {
517  //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
518  if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod))
519  {
520  fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
521  }
522  }// loop over super modules
523 
524  fEMCALGeoMatrixSet = kTRUE;//At least one, so good
525 
526  }//ESD as input
527  else
528  {
529  AliDebug(1,"Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file");
530  }//AOD as input
531  }//Get matrix from data
532  else if(gGeoManager)
533  {
534  fEMCALGeoMatrixSet = kTRUE;
535  }
536  }//EMCAL geo && no geoManager
537 
538  //Get the PHOS transformation geometry matrices from ESD
540  {
542  {
543  AliInfo("Load user defined PHOS geometry matrices");
544 
545  // OADB if available
546  AliOADBContainer geomContainer("phosGeo");
547  geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
548  TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");
549 
550  for(Int_t mod = 0 ; mod < 5 ; mod++)
551  {
552  if (!fPHOSMatrix[mod]) // Get it from OADB
553  {
554  AliDebug(1,Form("PHOS matrices module %d, %p",mod,((TGeoHMatrix*) matPHOS->At(mod))));
555  //((TGeoHMatrix*) matPHOS->At(mod))->Print();
556 
557  fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
558  }
559 
560  // Set it, if it exists
561  if(fPHOSMatrix[mod])
562  {
563  if(fDebug > 1 )
564  fPHOSMatrix[mod]->Print();
565 
566  fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
567  }
568  }// SM loop
569 
570  fPHOSGeoMatrixSet = kTRUE;//At least one, so good
571 
572  }//Load matrices
573  else if (!gGeoManager)
574  {
575  AliDebug(1,"Load PHOS misalignment matrices.");
576  if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
577  {
578  for(Int_t mod = 0; mod < 5; mod++)
579  {
580  if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod))
581  {
582  //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
583  fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
584  }
585  } // loop over modules
586 
587  fPHOSGeoMatrixSet = kTRUE; //At least one so good
588  } // ESD as input
589  else
590  {
591  AliDebug(1,"Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file");
592  } // AOD as input
593  } // get matrix from data
594  else if(gGeoManager)
595  {
596  fPHOSGeoMatrixSet = kTRUE;
597  }
598  }//PHOS geo and geoManager was not set
599 }
600 
601 //______________________________________________________________________________________
604 //______________________________________________________________________________________
606 {
607  Bool_t areNeighbours = kFALSE ;
608 
609  Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
610  Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
611 
612  Int_t rowdiff = 0, coldiff = 0;
613 
614  Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
615  Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
616 
617  if(calo==AliFiducialCut::kEMCAL && nSupMod1!=nSupMod2)
618  {
619  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
620  // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
621  if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
622  else icol2+=AliEMCALGeoParams::fgkEMCALCols;
623  }
624 
625  rowdiff = TMath::Abs( irow1 - irow2 ) ;
626  coldiff = TMath::Abs( icol1 - icol2 ) ;
627 
628  //if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
629  if ((coldiff + rowdiff == 1 ))
630  areNeighbours = kTRUE ;
631 
632  return areNeighbours;
633 }
634 
635 //_________________________________________________________________________________________
638 //_________________________________________________________________________________________
640  AliVCluster* cluster)
641 {
642  Int_t iSupMod = -1;
643  Int_t iSM0 = -1;
644  Int_t iTower = -1;
645  Int_t iIphi = -1;
646  Int_t iIeta = -1;
647  Int_t iphi = -1;
648  Int_t ieta = -1;
649 
650  for(Int_t iDigit = 0; iDigit < cluster->GetNCells(); iDigit++)
651  {
652  // Get from the absid the supermodule, tower and eta/phi numbers
653  geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
654  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
655 
656  // Check if there are cells of different SM
657  if (iDigit == 0 ) iSM0 = iSupMod;
658  else if(iSupMod != iSM0)
659  {
660  if(iSupMod > 11 && iSupMod < 18)
661  AliWarning(Form("Cluster shared in 2 DCal: SM%d, SM%d??",iSupMod,iSM0));
662 
663  return kTRUE;
664  }
665  }
666 
667  return kFALSE;
668 }
669 
670 //______________________________________________________________________________
673 //______________________________________________________________________________
675  AliVCaloCells* cells) const
676 {
677  //If the distance to the border is 0 or negative just exit accept all clusters
678 
679  if ( cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
680 
681  if ( cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
682 
683  Int_t absIdMax = -1;
684  Float_t ampMax = -1;
685 
686  for(Int_t i = 0; i < cluster->GetNCells() ; i++)
687  {
688  Int_t absId = cluster->GetCellAbsId(i) ;
689  Float_t amp = cells->GetCellAmplitude(absId);
690 
691  if(amp > ampMax)
692  {
693  ampMax = amp;
694  absIdMax = absId;
695  }
696  }
697 
698  AliDebug(1,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f",
699  absIdMax, ampMax, cluster->E()));
700 
701  if ( absIdMax == -1 ) return kFALSE;
702 
703  // Check if the cell is close to the borders:
704  Bool_t okrow = kFALSE;
705  Bool_t okcol = kFALSE;
706 
707  if ( cells->GetType() == AliVCaloCells::kEMCALCell )
708  {
709  // It should be the same as AliEMCALRecoUtils::CheckCellFiducialRegion()
710  // Why not calling it?
711 
712  Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
713 
714  fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
715 
716  fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
717 
718  if(iSM < 0 || iphi < 0 || ieta < 0 )
719  {
720  AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name",iSM,ieta,iphi));
721  }
722 
723  // Check rows/phi
724  Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
725 
726  if ( iSM < 10 || (iSM > 11 && iSM < 18) ) // Full EMCal (SM0-9) and DCal 2/3 (SM12-17)
727  {
728  if(iphi >= nborder && iphi < 24-nborder) okrow = kTRUE;
729  }
730  else // 1/3 SMs (SM10-11, SM18-19)
731  {
732  if(iphi >= nborder && iphi < 8-nborder) okrow = kTRUE;
733  }
734 
735  // Check columns/eta
736 
737  // Remove all borders if IsEMCALNoBorderAtEta0 or DCal SMs(12-17)
738  if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0() || (iSM > 11 && iSM < 18))
739  {
740  if(ieta > nborder && ieta < 48-nborder) okcol =kTRUE;
741  }
742  else // Do not remove borders close at eta = 0 for Full EMCal SMs and 1/3 EMCal
743  {
744  if ( iSM%2 == 0 )
745  {
746  if(ieta >= nborder) okcol = kTRUE;
747  }
748  else
749  {
750  if(ieta < 48-nborder) okcol = kTRUE;
751  }
752  } // eta 0 not checked
753 
754  AliDebug(1,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ? ok row %d, ok column %d",
755  nborder, ieta, iphi, iSM,okrow,okcol));
756 
757  }//EMCAL
758  else if ( cells->GetType() == AliVCaloCells::kPHOSCell )
759  {
760  Int_t relId[4];
761  Int_t irow = -1, icol = -1;
762  fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
763 
764  if (relId[1] != 0 ) return kFALSE; // skip CPV only PHOS
765 
766  irow = relId[2];
767  icol = relId[3];
768  //imod = relId[0]-1;
769 
770  if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
771  if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
772 
773  AliDebug(1,Form("PHOS Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ? ok row %d, ok column %d",
774  fNCellsFromPHOSBorder, icol, irow, relId[0]-1,okrow,okcol));
775  }//PHOS
776 
777  if (okcol && okrow) return kTRUE;
778  else return kFALSE;
779 }
780 
781 //________________________________________________________________________________________________________
784 //________________________________________________________________________________________________________
786 {
787  if (!fRemoveBadChannels) return kFALSE;
788 
789  //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
790  if(calorimeter == AliFiducialCut::kEMCAL && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
791  if(calorimeter == AliFiducialCut::kPHOS && !fPHOSBadChannelMap) return kFALSE;
792 
793  Int_t icol = -1;
794  Int_t irow = -1;
795  Int_t imod = -1;
796  for(Int_t iCell = 0; iCell<nCells; iCell++)
797  {
798  // Get the column and row
799  if ( calorimeter == AliFiducialCut::kEMCAL )
800  {
801  return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
802  }
803  else if ( calorimeter == AliFiducialCut::kPHOS )
804  {
805  Int_t relId[4];
806  fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
807 
808  if (relId[1] != 0 ) return kTRUE; // skip CPV only PHOS
809 
810  irow = relId[2];
811  icol = relId[3];
812  imod = relId[0]-1;
813 
814  if ( fPHOSBadChannelMap->GetEntries() <= imod ) continue;
815 
816  //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
817  if ( GetPHOSChannelStatus(imod, irow, icol) ) return kTRUE;
818  }
819  else return kFALSE;
820  } // cell cluster loop
821 
822  return kFALSE;
823 }
824 
825 //_______________________________________________________________
838 //_______________________________________________________________
840 (Bool_t bMC , Bool_t bExotic, Bool_t bNonLin,
841  Bool_t bRecalE, Bool_t bBad , Bool_t bRecalT, Int_t debug)
842 {
843  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - **** Start ***\n");
844 
845  // Exotic cells removal
846 
847  if(bExotic)
848  {
849  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Remove exotics in EMCAL\n");
850  fEMCALRecoUtils->SwitchOnRejectExoticCell() ;
851  fEMCALRecoUtils->SwitchOnRejectExoticCluster();
852 
853  // fEMCALRecoUtils->SetExoticCellDiffTimeCut(50); // If |t cell max - t cell in cross| > 50 do not add its energy, avoid
854  fEMCALRecoUtils->SetExoticCellFractionCut(0.97); // 1-Ecross/Ecell > 0.97 -> out
855  fEMCALRecoUtils->SetExoticCellMinAmplitudeCut(4.); // 4 GeV
856  }
857 
858  // Recalibration factors
859 
860  if(bRecalE && ! bMC)
861  {
862  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on energy recalibration in EMCAL\n");
863  fEMCALRecoUtils->SwitchOnRecalibration();
864  fEMCALRecoUtils->SwitchOnRunDepCorrection();
865  }
866 
867  // Remove EMCAL hot channels
868 
869  if(bBad)
870  {
871  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on bad channels removal in EMCAL\n");
872  fEMCALRecoUtils->SwitchOnBadChannelsRemoval();
873  fEMCALRecoUtils->SwitchOnDistToBadChannelRecalculation();
874  }
875 
876  // *** Time recalibration settings ***
877 
878  if(bRecalT && ! bMC)
879  {
880  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on time recalibration in EMCAL\n");
881  fEMCALRecoUtils->SwitchOnTimeRecalibration();
882  fEMCALRecoUtils->SwitchOnL1PhaseInTimeRecalibration() ;
883  }
884 
885  // Recalculate position with method
886 
887  fEMCALRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
888 
889  // Non linearity
890 
891  if( bNonLin )
892  {
893  if(!bMC)
894  {
895  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kBeamTestCorrected xxx\n");
896  fEMCALRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrectedv3);
897  }
898  else
899  {
900  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kPi0MCv3 xxx\n");
901  fEMCALRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MCv3);
902  }
903  }
904  else
905  {
906  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx DON'T SET Non linearity correction xxx\n");
907  fEMCALRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kNoCorrection);
908  }
909 
910 }
911 
912 //_______________________________________________________________
914 //_______________________________________________________________
916 {
917  clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
918 }
919 
920 //______________________________________________________________________________________
928 //______________________________________________________________________________________
929 Float_t AliCalorimeterUtils::GetECross(Int_t absID, AliVCaloCells* cells, Int_t bc)
930 {
931  if ( cells->IsEMCAL() )
932  {
933  Double_t tcell = cells->GetCellTime(absID);
934 
935  return fEMCALRecoUtils->GetECross(absID,tcell,cells,bc);
936  }
937  else // PHOS
938  {
939  Int_t icol = -1, irow = -1,iRCU = -1;
940  Int_t imod = GetModuleNumberCellIndexes(absID, AliFiducialCut::kPHOS, icol, irow, iRCU);
941 
942  Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
943 
944  Int_t relId1[] = { imod+1, 0, irow+1, icol };
945  Int_t relId2[] = { imod+1, 0, irow-1, icol };
946  Int_t relId3[] = { imod+1, 0, irow , icol+1 };
947  Int_t relId4[] = { imod+1, 0, irow , icol-1 };
948 
949  fPHOSGeo->RelToAbsNumbering(relId1, absId1);
950  fPHOSGeo->RelToAbsNumbering(relId2, absId2);
951  fPHOSGeo->RelToAbsNumbering(relId3, absId3);
952  fPHOSGeo->RelToAbsNumbering(relId4, absId4);
953 
954  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
955 
956  if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
957  if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
958  if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
959  if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
960 
961  return ecell1+ecell2+ecell3+ecell4;
962  }
963 }
964 
965 //_______________________________________________________________
976 //______________________________________________________________________________
977 void AliCalorimeterUtils::GetEMCALSubregion(AliVCluster * clus, AliVCaloCells * cells,
978  Int_t & regEta, Int_t & regPhi) const
979 {
980  regEta = regPhi = -1 ;
981 
982  if(!clus->IsEMCAL()) return ;
983 
984  Int_t icol = -1, irow = -1, iRCU = -1;
985  Float_t clusterFraction = 0;
986 
987  Int_t absId = GetMaxEnergyCell(cells,clus,clusterFraction);
988 
989  Int_t sm = GetModuleNumberCellIndexes(absId,AliFiducialCut::kEMCAL,icol,irow,iRCU);
990 
991  // Shift by 48 to for impair SM
992  if( sm%2 == 1) icol+=AliEMCALGeoParams::fgkEMCALCols;
993 
994  // Avoid borders
995  if(icol < 2 || icol > 93 || irow < 2 || irow > 21) return;
996 
997  //
998  // Eta regions
999  //
1000 
1001  // Region 0: center of SM ~0.18<|eta|<0.55
1002  if ( icol > 9 && icol < 34 ) regEta = 0;
1003  else if ( icol > 62 && icol < 87 ) regEta = 0;
1004 
1005  // Region 3: frames ~0.1<|eta|<~0.22 ~0.51<|eta|<0.62
1006 
1007  else if ( icol <= 9 && icol >= 5 ) regEta = 3;
1008  else if ( icol <= 38 && icol >= 34 ) regEta = 3;
1009  else if ( icol <= 62 && icol >= 58 ) regEta = 3;
1010  else if ( icol <= 91 && icol >= 87 ) regEta = 3;
1011 
1012  // Region 1: |eta| < ~0.15
1013 
1014  else if ( icol < 58 && icol > 38 ) regEta = 1 ;
1015 
1016  // Region 2: |eta| > ~0.6
1017 
1018  else regEta = 2 ;
1019 
1020  //
1021  // Phi regions
1022  //
1023 
1024  if ( irow >= 2 && irow <= 5 ) regPhi = 0; // External
1025  else if ( irow >= 18 && irow <= 21 ) regPhi = 0; // External
1026  else if ( irow >= 6 && irow <= 9 ) regPhi = 1; // Mid
1027  else if ( irow >= 14 && irow <= 17 ) regPhi = 1; // Mid
1028  else regPhi = 2; //10-13 Central
1029 
1030 }
1031 
1032 //________________________________________________________________________________________
1039 //________________________________________________________________________________________
1041 {
1042  // Get SM number
1043  Int_t sm = fEMCALGeo->GetSuperModuleNumber(absId);
1044 
1045  // Get first absId of SM
1046  Int_t absIdSMMin = fEMCALGeo->GetAbsCellIdFromCellIndexes(sm,0,1); // for absIds, first is in col 1, module/tower ordering map ...
1047 
1048  // Get reference n and correlated 3
1049  for(Int_t k = 0; k < 4; k++ )
1050  {
1051  for(Int_t p = 0; p < 72; p++ )
1052  {
1053  Int_t n = absIdSMMin + 2*k + 16 *p;
1054 
1055  if ( absId == n || absId == n+1 ||
1056  absId == n+8 || absId == n+9 )
1057  {
1058  absIdCorr[0] = n ;
1059  absIdCorr[1] = n+1 ;
1060  absIdCorr[2] = n+8 ;
1061  absIdCorr[3] = n+9 ;
1062 
1063  //printf("n=%d, n+1=%d, n+8=%d, n+9=%d\n",
1064  // absIdCorr[0],absIdCorr[1],absIdCorr[2],absIdCorr[3]);
1065 
1066  return kTRUE;
1067  }
1068  }
1069  }
1070 
1071  // Not found;
1072  absIdCorr[0] = -1 ;
1073  absIdCorr[1] = -1 ;
1074  absIdCorr[2] = -1 ;
1075  absIdCorr[3] = -1 ;
1076 
1077  return kFALSE;
1078 }
1079 
1080 //________________________________________________________________________________________
1089 //________________________________________________________________________________________
1091  Int_t & rowDiff, Int_t & colDiff) const
1092 {
1093  rowDiff = -100;
1094  colDiff = -100;
1095 
1096  if(absId1 == absId2) return kFALSE;
1097 
1098  // Check if in same SM, if not for sure not same TCard
1099  Int_t sm1 = fEMCALGeo->GetSuperModuleNumber(absId1);
1100  Int_t sm2 = fEMCALGeo->GetSuperModuleNumber(absId2);
1101  if ( sm1 != sm2 ) return kFALSE ;
1102 
1103  // Get the column and row of each absId
1104  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1105 
1106  Int_t col1, row1;
1107  fEMCALGeo->GetCellIndex(absId1,sm1,iTower,iIphi,iIeta);
1108  fEMCALGeo->GetCellPhiEtaIndexInSModule(sm1,iTower,iIphi, iIeta,row1,col1);
1109 
1110  Int_t col2, row2;
1111  fEMCALGeo->GetCellIndex(absId2,sm2,iTower,iIphi,iIeta);
1112  fEMCALGeo->GetCellPhiEtaIndexInSModule(sm2,iTower,iIphi, iIeta,row2,col2);
1113 
1114  Int_t row0 = Int_t(row1-row1%8);
1115  Int_t col0 = Int_t(col1-col1%2);
1116 
1117  Int_t rowDiff0 = row2-row0;
1118  Int_t colDiff0 = col2-col0;
1119 
1120  rowDiff = row1-row2;
1121  colDiff = col1-col2;
1122 
1123  // TCard is made by 2x8 towers
1124  if ( colDiff0 >=0 && colDiff0 < 2 && rowDiff0 >=0 && rowDiff0 < 8 )
1125  {
1126 
1127 // printf("\t absId (%d,%d), sm %d; col (%d,%d), colDiff %d; row (%d,%d),rowDiff %d\n",
1128 // absId1 , absId2, sm1,
1129 // col1, col2, colDiff,
1130 // row1, row2, rowDiff);
1131  return kTRUE ;
1132  }
1133  else
1134  return kFALSE;
1135 }
1136 
1137 
1138 //________________________________________________________________________________________
1140 //________________________________________________________________________________________
1142  AliVCluster * clu,
1143  Float_t & clusterFraction) const
1144 {
1145  if( !clu || !cells )
1146  {
1147  AliInfo("Cluster or cells pointer is null!");
1148  return -1;
1149  }
1150 
1151  Double_t eMax =-1.;
1152  Double_t eTot = 0.;
1153  Double_t eCell =-1.;
1154  Float_t fraction = 1.;
1155  Float_t recalFactor = 1.;
1156  Int_t cellAbsId =-1 , absId =-1 ;
1157  Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
1158 
1160  if(clu->IsPHOS()) calo = AliFiducialCut::kPHOS ;
1161 
1162  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1163  {
1164  cellAbsId = clu->GetCellAbsId(iDig);
1165 
1166  fraction = clu->GetCellAmplitudeFraction(iDig);
1167  if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1168 
1169  iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
1170 
1171  if(IsRecalibrationOn())
1172  {
1173  if(calo == AliFiducialCut::kEMCAL)
1174  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1175  else
1176  recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
1177  }
1178 
1179  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1180 
1181  if(eCell > eMax)
1182  {
1183  eMax = eCell;
1184  absId = cellAbsId;
1185  }
1186 
1187  eTot+=eCell;
1188 
1189  }// cell loop
1190 
1191  if(eTot > 0.1)
1192  clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
1193  else
1194  clusterFraction =1.;
1195 
1196  return absId;
1197 }
1198 
1199 //___________________________________________________________________________________
1203 //___________________________________________________________________________________
1204 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(AliVCluster* cluster,
1205  AliVEvent* event, Int_t index) const
1206 {
1207  AliVTrack *track = 0x0;
1208 
1209  //
1210  // EMCAL case only when matching is recalculated
1211  //
1212  if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
1213  {
1214  Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
1215  //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
1216 
1217  if(trackIndex < 0 )
1218  AliInfo(Form("Wrong track index %d, from recalculation", trackIndex));
1219  else
1220  track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
1221 
1222  return track ;
1223  }
1224 
1225  //
1226  // Normal case, get info from ESD or AOD
1227  //
1228 
1229  // No tracks matched
1230  if( cluster->GetNTracksMatched() < 1 ) return 0x0;
1231 
1232  // At least one match
1233  Int_t iTrack = 0; // only one match for AODs with index 0.
1234 
1235  // ESDs
1236  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
1237  {
1238  if( index >= 0 ) iTrack = index;
1239  else iTrack = ((AliESDCaloCluster*)cluster)->GetTracksMatched()->At(0); //cluster->GetTrackMatchedIndex();
1240 
1241  track = dynamic_cast<AliVTrack*> ( event->GetTrack(iTrack) );
1242  }
1243  else // AODs
1244  {
1245  if( index > 0 ) iTrack = index;
1246 
1247  track = dynamic_cast<AliVTrack*>( cluster->GetTrackMatched(iTrack) );
1248  }
1249 
1250  return track ;
1251 }
1252 
1258 {
1259  if( eCluster <= 0 || eCluster < eCell )
1260  {
1261  AliWarning(Form("Bad values eCell=%f, eCluster %f",eCell,eCluster));
1262  return 1;
1263  }
1264 
1265  Float_t frac = eCell / eCluster;
1266 
1267  Float_t correction = fMCECellClusFracCorrParam[0] +
1268  TMath::Exp( frac*fMCECellClusFracCorrParam[2]+fMCECellClusFracCorrParam[1] ) +
1269  fMCECellClusFracCorrParam[3]/TMath::Sqrt(frac);
1270 
1271 // printf("AliCalorimeterUtils::GetMCECellClusFracCorrection(eCell=%f, eCluster %f, frac %f) = %f\n",eCell, eCluster, frac, correction);
1272 // printf("\t %2.2f + TMath::Exp( %2.3f*%2.2f + %2.2f ) + %2.2f/TMath::Sqrt(%2.3f)) = %f\n",
1273 // fMCECellClusFracCorrParam[0],frac,fMCECellClusFracCorrParam[2],fMCECellClusFracCorrParam[1],fMCECellClusFracCorrParam[3], frac, correction);
1274 
1275  return correction;
1276 }
1277 
1278 //_____________________________________________________________________________________________________
1280 //_____________________________________________________________________________________________________
1281 Int_t AliCalorimeterUtils::GetModuleNumber(AliCaloTrackParticle * particle, AliVEvent * inputEvent) const
1282 {
1283  Int_t absId = -1;
1284 
1285  if(particle->GetDetectorTag()==AliFiducialCut::kEMCAL)
1286  {
1287  fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
1288 
1289  AliDebug(2,Form("EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d",
1290  particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId)));
1291 
1292  return fEMCALGeo->GetSuperModuleNumber(absId) ;
1293  } // EMCAL
1294  else if ( particle->GetDetectorTag() == AliFiducialCut::kPHOS )
1295  {
1296  // In case we use the MC reader, the input are TParticles,
1297  // in this case use the corresponing method in PHOS Geometry to get the particle.
1298  if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
1299  {
1300  Int_t mod =-1;
1301  Double_t z = 0., x=0.;
1302 
1303  TParticle* primary = ((AliMCEvent*)inputEvent)->Particle(particle->GetCaloLabel(0));
1304 
1305  if(primary)
1306  {
1307  fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
1308  }
1309  else
1310  {
1311  AliFatal("Primary not available, stop!");
1312  }
1313  return mod;
1314  }
1315  // Input are ESDs or AODs, get the PHOS module number like this.
1316  else
1317  {
1318  //FIXME
1319  //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
1320  //return GetModuleNumber(cluster);
1321  //MEFIX
1322  return -1;
1323  }
1324  } // PHOS
1325 
1326  return -1;
1327 }
1328 
1329 //_____________________________________________________________________
1331 //_____________________________________________________________________
1332 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
1333 {
1334  if(!cluster)
1335  {
1336  AliDebug(1,"AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
1337 
1338  return -1;
1339  }
1340 
1341  if ( cluster->GetNCells() <= 0 ) return -1;
1342 
1343  Int_t absId = cluster->GetCellAbsId(0);
1344 
1345  if ( absId < 0 ) return -1;
1346 
1347  if( cluster->IsEMCAL() )
1348  {
1349  AliDebug(2,Form("EMCAL absid %d, SuperModule %d",absId, fEMCALGeo->GetSuperModuleNumber(absId)));
1350 
1351  return fEMCALGeo->GetSuperModuleNumber(absId) ;
1352  } // EMCAL
1353  else if ( cluster->IsPHOS() )
1354  {
1355  Int_t relId[4];
1356  fPHOSGeo->AbsToRelNumbering(absId,relId);
1357 
1358  if (relId[1] != 0 ) return -1; // skip CPV only PHOS
1359 
1360  AliDebug(2,Form("PHOS absid %d Module %d",absId, relId[0]-1));
1361 
1362  return relId[0]-1;
1363  } // PHOS
1364 
1365  return -1;
1366 }
1367 
1368 //___________________________________________________________________________________________________
1370 //___________________________________________________________________________________________________
1372  Int_t & icol, Int_t & irow, Int_t & iRCU) const
1373 {
1374  Int_t imod = -1;
1375 
1376  if ( absId < 0 ) return -1 ;
1377 
1378  if ( calo == AliFiducialCut::kEMCAL || calo == AliFiducialCut::kDCAL)
1379  {
1380  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1381  fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1382  fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1383 
1384  if(imod < 0 || irow < 0 || icol < 0 )
1385  {
1386  AliFatal(Form("Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name",imod,icol,irow));
1387  }
1388 
1389  // In case of DCal C side, shift columns to match offline/online numbering
1390  // when calculating the DDL for Run2
1391  // See AliEMCALRawUtils::Digits2Raw and Raw2Digits.
1392  Int_t ico2 = icol ;
1393  if ( imod == 13 || imod == 15 || imod == 17 ) ico2 += 16;
1394 
1395  // RCU / DDL
1396  if(imod < 10 || (imod > 11 && imod < 18)) // (EMCAL Full || DCAL 2/3)
1397  {
1398  // RCU0 / DDL0
1399  if ( 0 <= irow && irow < 8 ) iRCU = 0; // first cable row
1400  else if ( 8 <= irow && irow < 16 &&
1401  0 <= ico2 && ico2 < 24 ) iRCU = 0; // first half;
1402  //second cable row
1403 
1404  // RCU1 / DDL1
1405  else if ( 8 <= irow && irow < 16 &&
1406  24 <= ico2 && ico2 < 48 ) iRCU = 1; // second half;
1407  //second cable row
1408  else if ( 16 <= irow && irow < 24 ) iRCU = 1; // third cable row
1409 
1410  if ( imod%2 == 1 ) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1411  }
1412  else
1413  {
1414  // 1/3 SM have one single SRU, just assign RCU/DDL 0
1415  iRCU = 0 ;
1416  }
1417 
1418  if ( iRCU < 0 )
1419  AliFatal(Form("Wrong EMCAL RCU number = %d", iRCU));
1420 
1421  return imod ;
1422  } // EMCAL
1423  else // PHOS
1424  {
1425  Int_t relId[4];
1426  fPHOSGeo->AbsToRelNumbering(absId,relId);
1427 
1428  if (relId[1] != 0 ) return -1; // skip CPV only PHOS
1429 
1430  irow = relId[2];
1431  icol = relId[3];
1432  imod = relId[0]-1;
1433  iRCU= (Int_t)(relId[2]-1)/16 ;
1434 
1435  //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
1436 
1437  if ( iRCU >= 4 )
1438  AliFatal(Form("Wrong PHOS RCU number = %d", iRCU));
1439 
1440  return imod;
1441  } // PHOS
1442 
1443  return -1;
1444 }
1445 
1446 //___________________________________________________________________________________________________
1449 //___________________________________________________________________________________________________
1451  Int_t & icol , Int_t & irow , Int_t & iRCU,
1452  Int_t & icolAbs, Int_t & irowAbs) const
1453 {
1454  Int_t imod = GetModuleNumberCellIndexes(absId, calo, icol, irow,iRCU);
1455 
1456  icolAbs = icol;
1457  irowAbs = irow;
1458 
1459  //
1460  // PHOS
1461  //
1462  if(calo == AliFiducialCut::kPHOS)
1463  {
1464  irowAbs = irow + 64 * imod;
1465 
1466  return imod;
1467  }
1468  //
1469  // EMCal/DCal
1470  //
1471  else
1472  {
1473  //
1474  // Shift collumns in even SM
1475  Int_t shiftEta = 48;
1476 
1477  // Shift collumn even more due to smaller acceptance of DCal collumns
1478  if ( imod > 11 && imod < 18) shiftEta+=48/3;
1479 
1480  icolAbs = (imod % 2) ? icol + shiftEta : icol;
1481 
1482  //
1483  // Shift rows per sector
1484  irowAbs = irow + 24 * Int_t(imod / 2);
1485 
1486  // Shift row less due to smaller acceptance of SM 10 and 11 to count DCal rows
1487  if ( imod > 11 && imod < 20) irowAbs -= (2*24 / 3);
1488 
1489  return imod ;
1490  }
1491 }
1492 
1493 
1494 //___________________________________________________________________________________________
1496 //___________________________________________________________________________________________
1497 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1498 {
1499  const Int_t nc = cluster->GetNCells();
1500 
1501  Int_t absIdList[nc];
1502  Float_t maxEList[nc];
1503 
1504  Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1505 
1506  return nMax;
1507 }
1508 
1509 //___________________________________________________________________________________________
1511 //___________________________________________________________________________________________
1512 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1513  Int_t *absIdList, Float_t *maxEList)
1514 {
1515  Int_t iDigitN = 0 ;
1516  Int_t iDigit = 0 ;
1517  Int_t absId1 = -1 ;
1518  Int_t absId2 = -1 ;
1519  const Int_t nCells = cluster->GetNCells();
1520 
1521  Float_t eCluster = RecalibrateClusterEnergy(cluster, cells);// recalculate cluster energy, avoid non lin correction.
1522 
1523  Float_t simuTotWeight = 0;
1525  {
1526  simuTotWeight = RecalibrateClusterEnergyWeightCell(cluster, cells,eCluster);// same but apply a weight
1527  simuTotWeight/= eCluster;
1528  }
1529 
1531  if(!cluster->IsEMCAL()) calorimeter = AliFiducialCut::kPHOS;
1532 
1533  //printf("cluster : ncells %d \n",nCells);
1534 
1535  Float_t emax = 0;
1536  Int_t idmax =-1;
1537  for(iDigit = 0; iDigit < nCells ; iDigit++)
1538  {
1539  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
1540  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1541  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1542 
1544  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1545 
1546  if( en > emax )
1547  {
1548  emax = en ;
1549  idmax = absIdList[iDigit] ;
1550  }
1551  //Int_t icol = -1, irow = -1, iRCU = -1;
1552  //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1553  //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
1554  }
1555 
1556  for(iDigit = 0 ; iDigit < nCells; iDigit++)
1557  {
1558  if( absIdList[iDigit] >= 0 )
1559  {
1560  absId1 = cluster->GetCellsAbsId()[iDigit];
1561 
1562  Float_t en1 = cells->GetCellAmplitude(absId1);
1563  RecalibrateCellAmplitude(en1,calorimeter,absId1);
1564 
1566  en1*=GetMCECellClusFracCorrection(en1,eCluster)/simuTotWeight;
1567 
1568  //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
1569 
1570  for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1571  {
1572  absId2 = cluster->GetCellsAbsId()[iDigitN] ;
1573 
1574  if(absId2==-1 || absId2==absId1) continue;
1575 
1576  //printf("\t %d : absIDj %d\n",iDigitN, absId2);
1577 
1578  Float_t en2 = cells->GetCellAmplitude(absId2);
1579  RecalibrateCellAmplitude(en2,calorimeter,absId2);
1580 
1582  en2*=GetMCECellClusFracCorrection(en2,eCluster)/simuTotWeight;
1583 
1584  //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
1585 
1586  if ( AreNeighbours(calorimeter, absId1, absId2) )
1587  {
1588  // printf("\t \t Neighbours \n");
1589  if ( en1 > en2 )
1590  {
1591  absIdList[iDigitN] = -1 ;
1592  //printf("\t \t indexN %d not local max\n",iDigitN);
1593  // but may be digit too is not local max ?
1594  if(en1 < en2 + fLocMaxCutEDiff) {
1595  //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
1596  absIdList[iDigit] = -1 ;
1597  }
1598  }
1599  else
1600  {
1601  absIdList[iDigit] = -1 ;
1602  //printf("\t \t index %d not local max\n",iDigitN);
1603  // but may be digitN too is not local max ?
1604  if(en1 > en2 - fLocMaxCutEDiff)
1605  {
1606  absIdList[iDigitN] = -1 ;
1607  //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
1608  }
1609  }
1610  } // if Are neighbours
1611  //else printf("\t \t NOT Neighbours \n");
1612  } // while digitN
1613  } // slot not empty
1614  } // while digit
1615 
1616  iDigitN = 0 ;
1617  for(iDigit = 0; iDigit < nCells; iDigit++)
1618  {
1619  if( absIdList[iDigit] >= 0 )
1620  {
1621  absIdList[iDigitN] = absIdList[iDigit] ;
1622 
1623  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1624  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1625 
1627  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1628 
1629  if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1630 
1631  maxEList[iDigitN] = en ;
1632 
1633  //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1634  iDigitN++ ;
1635  }
1636  }
1637 
1638  if ( iDigitN == 0 )
1639  {
1640  AliDebug(1,Form("No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f",
1641  idmax,emax,cluster->E()));
1642  iDigitN = 1 ;
1643  maxEList[0] = emax ;
1644  absIdList[0] = idmax ;
1645  }
1646 
1647 
1648  AliDebug(1,Form("In cluster E %2.2f (wth non lin. %2.2f), M02 %2.2f, M20 %2.2f, N maxima %d",
1649  cluster->E(),eCluster, cluster->GetM02(),cluster->GetM20(), iDigitN));
1650 
1651 // if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1652 // {
1653 // printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1654 // }
1655 
1656  return iDigitN ;
1657 }
1658 
1659 //____________________________________
1661 //____________________________________
1663 {
1664  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1665  {
1666  AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null");
1667  return TString("");
1668  }
1669 
1670  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1671  {
1672  AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null");
1673  return TString("");
1674  }
1675 
1676  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1677  if (pass.Contains("ass1")) return TString("pass1");
1678  else if (pass.Contains("ass2")) return TString("pass2");
1679  else if (pass.Contains("ass3")) return TString("pass3");
1680  else if (pass.Contains("ass4")) return TString("pass4");
1681  else if (pass.Contains("ass5")) return TString("pass5");
1682  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1683  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1684  {
1685  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1686  return TString("pass1");
1687  }
1688  else if (pass.Contains("LHC14a1a"))
1689  {
1690  AliInfo("Check that Energy calibration was enabled for this MC production in the tender, clusterizer or here!!");
1691 
1692  return TString("LHC14a1a");
1693  }
1694 
1695  // No condition fullfilled, give a default value
1696  AliInfo("Pass number string not found");
1697  return TString("");
1698 }
1699 
1700 //________________________________________
1702 //________________________________________
1704 {
1705  fEMCALGeoName = "";
1706  fPHOSGeoName = "";
1707 
1708  fEMCALGeoMatrixSet = kFALSE;
1709  fPHOSGeoMatrixSet = kFALSE;
1710 
1711  fRemoveBadChannels = kFALSE;
1712 
1714 
1715  fLocMaxCutE = 0.1 ;
1716  fLocMaxCutEDiff = 0.0 ;
1717 
1718  // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1719  // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1720  // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1721  // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1722  // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1723  // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
1724 
1725  fOADBSet = kFALSE;
1726  fOADBForEMCAL = kTRUE ;
1727  fOADBForPHOS = kFALSE;
1728 
1729  fOADBFilePathEMCAL = "" ; // $ALICE_PHYSICS/OADB/EMCAL
1730  fOADBFilePathPHOS = "$ALICE_PHYSICS/OADB/PHOS" ;
1731 
1732  fImportGeometryFromFile = kTRUE;
1734 
1735  fNSuperModulesUsed = 22;
1736 
1737  fMCECellClusFracCorrParam[0] = 0.78;
1738  fMCECellClusFracCorrParam[1] =-1.8;
1739  fMCECellClusFracCorrParam[2] =-6.3;
1740  fMCECellClusFracCorrParam[3] = 0.014;
1741 }
1742 
1743 //_____________________________________________________
1745 //_____________________________________________________
1747 {
1748  AliDebug(1,"Init bad channel map");
1749 
1750  // In order to avoid rewriting the same histograms
1751  Bool_t oldStatus = TH1::AddDirectoryStatus();
1752  TH1::AddDirectory(kFALSE);
1753 
1754  fPHOSBadChannelMap = new TObjArray(5);
1755  for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),
1756  Form("PHOS_BadMap_mod%d",i),
1757  64, 0, 64, 56, 0, 56));
1758 
1759  fPHOSBadChannelMap->SetOwner(kTRUE);
1760  fPHOSBadChannelMap->Compress();
1761 
1762  //In order to avoid rewriting the same histograms
1763  TH1::AddDirectory(oldStatus);
1764 }
1765 
1766 //______________________________________________________
1768 //______________________________________________________
1770 {
1771  AliDebug(1,"Init recalibration map");
1772 
1773  // In order to avoid rewriting the same histograms
1774  Bool_t oldStatus = TH1::AddDirectoryStatus();
1775  TH1::AddDirectory(kFALSE);
1776 
1778  for (int i = 0; i < 5; i++)
1779  {
1780  fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),
1781  Form("PHOSRecalFactors_Mod%d",i),
1782  64, 0, 64, 56, 0, 56));
1783  }
1784 
1785  // Init the histograms with 1
1786  for (Int_t m = 0; m < 5; m++)
1787  {
1788  for (Int_t i = 0; i < 56; i++)
1789  {
1790  for (Int_t j = 0; j < 64; j++)
1791  {
1793  }
1794  }
1795  }
1796 
1797  fPHOSRecalibrationFactors->SetOwner(kTRUE);
1798  fPHOSRecalibrationFactors->Compress();
1799 
1800  // In order to avoid rewriting the same histograms
1801  TH1::AddDirectory(oldStatus);
1802 }
1803 
1808 {
1809  if (fEMCALGeo) return;
1810 
1811  AliDebug(1,Form(" for run=%d",fRunNumber));
1812 
1813  if(fEMCALGeoName=="")
1814  {
1815  fEMCALGeo = AliEMCALGeometry::GetInstanceFromRunNumber(fRunNumber);
1816  AliInfo(Form("Get EMCAL geometry name to <%s> for run %d",fEMCALGeo->GetName(),fRunNumber));
1817  }
1818  else
1819  {
1820  fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
1821  AliInfo(Form("Set EMCAL geometry name to <%s>",fEMCALGeoName.Data()));
1822  }
1823 
1824  // Init geometry, I do not like much to do it like this ...
1825  if(fImportGeometryFromFile && !gGeoManager)
1826  {
1827  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1828  {
1829  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1830  if (fRunNumber < 140000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2010.root").data();
1831  else if(fRunNumber < 171000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2011.root").data();
1832  else if(fRunNumber < 198000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2012.root").data(); // 2012-2013
1833  else fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2015.root").data(); // >= 2015
1834  }
1835 
1836  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1837 
1838  TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1839  }
1840  else if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1841 }
1842 
1847 {
1848  if (fPHOSGeo) return;
1849 
1850  AliDebug(1,Form(" for run=%d",fRunNumber));
1851 
1852  if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1853 
1854  fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
1855 
1856  //if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1857 }
1858 
1859 //______________________________________________________________________________________________
1860 // Check that a MC ESD is in the calorimeter acceptance
1861 //______________________________________________________________________________________________
1863 {
1864  if(!particle || (calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS)) return kFALSE ;
1865 
1866  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
1868  {
1869  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1870  return kFALSE ;
1871  }
1872 
1873  if(calo == AliFiducialCut::kPHOS )
1874  {
1875  Int_t mod = 0 ;
1876  Double_t x = 0, z = 0 ;
1877  return GetPHOSGeometry()->ImpactOnEmc( particle, mod, z, x);
1878  }
1879  else if(calo == AliFiducialCut::kEMCAL)
1880  {
1881  Int_t absID = 0 ;
1882  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(),absID);
1883  if(ok)
1884  {
1885  Int_t icol = -1, irow = -1, iRCU = -1;
1886  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1887  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1888  if(status > 0) ok = kFALSE;
1889  }
1890 
1891  return ok ;
1892  }
1893 
1894  return kFALSE ;
1895 }
1896 
1897 //______________________________________________________________________________________________________
1899 //______________________________________________________________________________________________________
1901 {
1902  if(!particle || (calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS)) return kFALSE ;
1903 
1904  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
1906  {
1907  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1908  return kFALSE ;
1909  }
1910 
1911  Float_t phi = particle->Phi();
1912  if(phi < 0) phi+=TMath::TwoPi();
1913 
1914  if(calo == AliFiducialCut::kPHOS )
1915  {
1916  Int_t mod = 0 ;
1917  Double_t x = 0, z = 0 ;
1918  Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1919  return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), phi, mod, z, x) ;
1920  }
1921  else if(calo == AliFiducialCut::kEMCAL)
1922  {
1923  Int_t absID = 0 ;
1924  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),phi,absID);
1925  if(ok)
1926  {
1927  Int_t icol = -1, irow = -1, iRCU = -1;
1928  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1929  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1930  if(status > 0) ok = kFALSE;
1931  }
1932 
1933  return ok ;
1934  }
1935 
1936  return kFALSE ;
1937 }
1938 
1939 //______________________________________________________________________________________________________
1941 //______________________________________________________________________________________________________
1943 {
1944  if(!particle || (calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS)) return kFALSE ;
1945 
1946  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
1948  {
1949  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1950  return kFALSE ;
1951  }
1952 
1953  Float_t phi = particle->Phi();
1954  if(phi < 0) phi+=TMath::TwoPi();
1955 
1956  if(calo == AliFiducialCut::kPHOS )
1957  {
1958  Int_t mod = 0 ;
1959  Double_t x = 0, z = 0 ;
1960  Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1961  return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), phi, mod, z, x) ;
1962  }
1963  else if(calo == AliFiducialCut::kEMCAL)
1964  {
1965  Int_t absID = 0 ;
1966  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),phi,absID);
1967  if(ok)
1968  {
1969  Int_t icol = -1, irow = -1, iRCU = -1;
1970  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1971  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1972  if(status > 0) ok = kFALSE;
1973  }
1974 
1975  return ok ;
1976  }
1977 
1978  return kFALSE ;
1979 }
1980 
1981 //_____________________________________________________________________________________________________
1982 // Check that a TLorentzVector is in the calorimeter acceptance, give the cell number where it hit.
1983 //_____________________________________________________________________________________________________
1985  Float_t phiOrg, Int_t & absID)
1986 {
1987  if(calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS) return kFALSE ;
1988 
1989  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
1991  {
1992  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1993  return kFALSE ;
1994  }
1995 
1996  Float_t phi = phiOrg;
1997  if(phi < 0) phi+=TMath::TwoPi();
1998 
1999  if(calo == AliFiducialCut::kPHOS )
2000  {
2001  Int_t mod = 0 ;
2002  Double_t x = 0, z = 0 ;
2003  Double_t vtx[]={0,0,0} ;
2004  return GetPHOSGeometry()->ImpactOnEmc(vtx, theta, phi, mod, z, x) ;
2005  }
2006  else if(calo == AliFiducialCut::kEMCAL)
2007  {
2008  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta,phi,absID);
2009  if(ok)
2010  {
2011  Int_t icol = -1, irow = -1, iRCU = -1;
2012  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
2013  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
2014  if(status > 0) ok = kFALSE;
2015  }
2016 
2017  return ok ;
2018  }
2019 
2020  return kFALSE ;
2021 }
2022 
2023 //_______________________________________________________________________
2026 //_______________________________________________________________________
2028 {
2029  Int_t icol = ieta;
2030  if ( iSM%2 ) icol+=48; // Impair SM, shift index [0-47] to [48-96]
2031 
2033  {
2034  for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
2035  {
2036  if(icol==fMaskCellColumns[imask]) return kTRUE;
2037  }
2038  }
2039 
2040  return kFALSE;
2041 }
2042 
2043 //_________________________________________________________
2045 //_________________________________________________________
2046 void AliCalorimeterUtils::Print(const Option_t * opt) const
2047 {
2048  if(! opt)
2049  return;
2050 
2051  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
2052  printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
2053  printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
2054  fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
2055  if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
2056  printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
2057  printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
2058  printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
2059  printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
2060 
2061  printf("Recalibrate time? %d, With L1 phase run by run? %d\n",IsTimeRecalibrationOn(),IsL1PhaseInTimeRecalibrationOn());
2062 
2063  printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
2064  printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
2065 
2066  printf(" \n") ;
2067 }
2068 
2069 //_____________________________________________________________________________________________
2071 //_____________________________________________________________________________________________
2073 {
2074  Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
2075  Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
2076 
2077  if (IsRecalibrationOn())
2078  {
2079  if(calo == AliFiducialCut::kPHOS)
2080  {
2081  amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
2082  }
2083  else
2084  {
2085  amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
2086  }
2087  }
2088 }
2089 
2090 //____________________________________________________________________________________________________
2092 //____________________________________________________________________________________________________
2094 {
2096  {
2097  GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
2098  }
2099 }
2100 
2101 
2102 //____________________________________________________________________________________________________
2104 //____________________________________________________________________________________________________
2105 void AliCalorimeterUtils::RecalibrateCellTimeL1Phase(Double_t & time, Int_t calo, Int_t iSM, Int_t bunchCrossNumber) const
2106 {
2108  {
2109  GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(iSM, bunchCrossNumber, time);
2110  }
2111 }
2112 
2113 
2114 //__________________________________________________________________________
2116 //__________________________________________________________________________
2118  AliVCaloCells * cells)
2119 {
2120  if ( !IsRecalibrationOn() ) return cluster->E();
2121 
2122  // Initialize some used variables
2123  Float_t frac = 0., energy = 0.;
2124 
2125  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
2126  UShort_t * index = cluster->GetCellsAbsId() ;
2127  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
2128  Int_t ncells = cluster->GetNCells();
2129 
2131  if(cluster->IsPHOS()) calo = AliFiducialCut::kPHOS ;
2132 
2133  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
2134  for(Int_t icell = 0; icell < ncells; icell++)
2135  {
2136  Int_t absId = index[icell];
2137 
2138  frac = fraction[icell];
2139  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
2140 
2141  Float_t amp = cells->GetCellAmplitude(absId);
2142  RecalibrateCellAmplitude(amp,calo, absId);
2143 
2144  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy: before cal %f; after cal %f",
2145  calo,frac,cells->GetCellAmplitude(absId),amp));
2146 
2147  energy += amp*frac;
2148  }
2149 
2150  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
2151 
2152  return energy;
2153 }
2154 
2155 //_______________________________________________________________________________________________________
2158 //_______________________________________________________________________________________________________
2160  AliVCaloCells * cells, Float_t energyOrg)
2161 {
2162  // Initialize some used variables
2163  Float_t frac = 0., energy = 0.;
2164 
2165  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
2166  UShort_t * index = cluster->GetCellsAbsId() ;
2167  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
2168  Int_t ncells = cluster->GetNCells();
2169 
2171  if(cluster->IsPHOS()) calo = AliFiducialCut::kPHOS ;
2172 
2173  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
2174  for(Int_t icell = 0; icell < ncells; icell++)
2175  {
2176  Int_t absId = index[icell];
2177 
2178  frac = fraction[icell];
2179  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
2180 
2181  Float_t amp = cells->GetCellAmplitude(absId);
2182  if ( IsRecalibrationOn() ) RecalibrateCellAmplitude(amp,calo, absId);
2183 
2184  amp*=GetMCECellClusFracCorrection(amp,energyOrg);
2185 
2186  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy %f",
2187  calo,frac,cells->GetCellAmplitude(absId)));
2188 
2189  energy += amp*frac;
2190  }
2191 
2192  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
2193 
2194  return energy;
2195 }
2196 
2197 //__________________________________________________________________________________________
2200 //__________________________________________________________________________________________
2201 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
2202 {
2203  fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
2204 }
2205 
2206 //________________________________________________________________________________
2213 //________________________________________________________________________________
2215  TObjArray* clusterArray,
2216  AliMCEvent* mc)
2217 {
2218  if (!fRecalculateMatching) return ;
2219 
2220  fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo,mc) ;
2221 
2222  Float_t dZ = 2000;
2223  Float_t dR = 2000;
2224 
2225  Int_t nClusters = event->GetNumberOfCaloClusters();
2226  if(clusterArray) nClusters = clusterArray->GetEntriesFast();
2227 
2228  AliVCluster * clus = 0;
2229 
2230  for (Int_t iclus = 0; iclus < nClusters ; iclus++)
2231  {
2232  if ( clusterArray ) clus = (AliVCluster*) clusterArray->At(iclus) ;
2233  else clus = event->GetCaloCluster(iclus) ;
2234 
2235  if (!clus->IsEMCAL()) continue ;
2236 
2237  //
2238  // Put track residuals in cluster
2239  //
2240  fEMCALRecoUtils->GetMatchedResiduals(clus->GetID(),dZ,dR);
2241 
2242  if ( TMath::Abs(clus->GetTrackDx()) < 500 )
2243  AliDebug(2,Form("Residuals (Old, New): z (%2.4f,%2.4f), x (%2.4f,%2.4f)\n",
2244  clus->GetTrackDz(),dZ,clus->GetTrackDx(),dR));
2245 
2246  clus->SetTrackDistance(dR,dZ);
2247 
2248  //
2249  // Remove old matches in cluster
2250  //
2251  if(clus->GetNTracksMatched() > 0)
2252  {
2253  if(!strcmp("AliESDCaloCluster",Form("%s",clus->ClassName())))
2254  {
2255  TArrayI arrayTrackMatched(0);
2256  ((AliESDCaloCluster*)clus)->AddTracksMatched(arrayTrackMatched);
2257  }
2258  else
2259  {
2260  for(Int_t iTrack = 0; iTrack < clus->GetNTracksMatched(); iTrack++)
2261  {
2262  ((AliAODCaloCluster*)clus)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)clus)->GetTrackMatched(iTrack));
2263  }
2264  }
2265  }
2266 
2267  //
2268  // Now put first track index in cluster.
2269  //
2270  Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(iclus);
2271  if ( trackIndex >= 0 )
2272  {
2273  if(!strcmp("AliESDCaloCluster",Form("%s",clus->ClassName())))
2274  {
2275  TArrayI arrayTrackMatched(1);
2276  arrayTrackMatched[0] = trackIndex;
2277  ((AliESDCaloCluster*)clus)->AddTracksMatched(arrayTrackMatched);
2278  }
2279  else
2280  {
2281  ((AliAODCaloCluster*)clus)->AddTrackMatched((TObject*)event->GetTrack(trackIndex));
2282  }
2283  }
2284 
2285  } // cluster loop
2286 }
2287 
2288 //___________________________________________________________________________
2301 //___________________________________________________________________________
2303  AliVCluster* cluster,
2304  AliVCaloCells* cells,
2305  //Float_t & e1, Float_t & e2,
2306  AliAODCaloCluster* cluster1,
2307  AliAODCaloCluster* cluster2,
2308  Int_t nMax, Int_t eventNumber)
2309 {
2310  TH2F* hClusterMap = 0 ;
2311  TH2F* hClusterLocMax = 0 ;
2312  TH2F* hCluster1 = 0 ;
2313  TH2F* hCluster2 = 0 ;
2314 
2315  if(fPlotCluster)
2316  {
2317  hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
2318  hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
2319  hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
2320  hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
2321  hClusterMap ->SetXTitle("column");
2322  hClusterMap ->SetYTitle("row");
2323  hClusterLocMax ->SetXTitle("column");
2324  hClusterLocMax ->SetYTitle("row");
2325  hCluster1 ->SetXTitle("column");
2326  hCluster1 ->SetYTitle("row");
2327  hCluster2 ->SetXTitle("column");
2328  hCluster2 ->SetYTitle("row");
2329  }
2330 
2332  if(cluster->IsPHOS())
2333  {
2334  calorimeter = AliFiducialCut::kPHOS;
2335  AliWarning("Not supported for PHOS yet");
2336  return;
2337  }
2338 
2339  const Int_t ncells = cluster->GetNCells();
2340  Int_t absIdList[ncells];
2341 
2342  Float_t e1 = 0, e2 = 0 ;
2343  Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
2344  Float_t eCluster = 0;
2345  Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
2346  for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
2347  {
2348  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
2349 
2350  Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
2351  RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
2352  eCluster+=ec;
2353 
2354  if(fPlotCluster)
2355  {
2356  //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
2357  sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
2358  if(sm > -1 && sm < 12) // just to avoid compilation warning
2359  {
2360  if(icol > maxCol) maxCol = icol;
2361  if(icol < minCol) minCol = icol;
2362  if(irow > maxRow) maxRow = irow;
2363  if(irow < minRow) minRow = irow;
2364  hClusterMap->Fill(icol,irow,ec);
2365  }
2366  }
2367  }
2368 
2369  // Init counters and variables
2370  Int_t ncells1 = 1 ;
2371  UShort_t absIdList1[9] ;
2372  Double_t fracList1 [9] ;
2373  absIdList1[0] = absId1 ;
2374  fracList1 [0] = 1. ;
2375 
2376  Float_t ecell1 = cells->GetCellAmplitude(absId1);
2377  RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
2378  e1 = ecell1;
2379 
2380  Int_t ncells2 = 1 ;
2381  UShort_t absIdList2[9] ;
2382  Double_t fracList2 [9] ;
2383  absIdList2[0] = absId2 ;
2384  fracList2 [0] = 1. ;
2385 
2386  Float_t ecell2 = cells->GetCellAmplitude(absId2);
2387  RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
2388  e2 = ecell2;
2389 
2390  if(fPlotCluster)
2391  {
2392  Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
2393  sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
2394  hClusterLocMax->Fill(icol1,irow1,ecell1);
2395  sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
2396  hClusterLocMax->Fill(icol2,irow2,ecell2);
2397  }
2398 
2399  // Very rough way to share the cluster energy
2400  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
2401  Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
2402  Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
2403 
2404  for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
2405  {
2406  Int_t absId = absIdList[iDigit];
2407 
2408  if ( absId==absId1 || absId==absId2 || absId < 0 ) continue;
2409 
2410  Float_t ecell = cells->GetCellAmplitude(absId);
2411  RecalibrateCellAmplitude(ecell, calorimeter, absId);
2412 
2413  if(AreNeighbours(calorimeter, absId1,absId ))
2414  {
2415  absIdList1[ncells1]= absId;
2416 
2417  if(AreNeighbours(calorimeter, absId2,absId ))
2418  {
2419  fracList1[ncells1] = shareFraction1;
2420  e1 += ecell*shareFraction1;
2421  }
2422  else
2423  {
2424  fracList1[ncells1] = 1.;
2425  e1 += ecell;
2426  }
2427 
2428  ncells1++;
2429 
2430  } // neigbour to cell1
2431 
2432  if(AreNeighbours(calorimeter, absId2,absId ))
2433  {
2434  absIdList2[ncells2]= absId;
2435 
2436  if(AreNeighbours(calorimeter, absId1,absId ))
2437  {
2438  fracList2[ncells2] = shareFraction2;
2439  e2 += ecell*shareFraction2;
2440  }
2441  else
2442  {
2443  fracList2[ncells2] = 1.;
2444  e2 += ecell;
2445  }
2446 
2447  ncells2++;
2448 
2449  } // neigbour to cell2
2450  }
2451 
2452  AliDebug(1,Form("N Local Max %d, Cluster energy = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2 %f, sum f12 = %f",
2453  nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2));
2454 
2455  cluster1->SetE(e1);
2456  cluster2->SetE(e2);
2457 
2458  cluster1->SetNCells(ncells1);
2459  cluster2->SetNCells(ncells2);
2460 
2461  cluster1->SetCellsAbsId(absIdList1);
2462  cluster2->SetCellsAbsId(absIdList2);
2463 
2464  cluster1->SetCellsAmplitudeFraction(fracList1);
2465  cluster2->SetCellsAmplitudeFraction(fracList2);
2466 
2467  // Correct linearity
2468  CorrectClusterEnergy(cluster1) ;
2469  CorrectClusterEnergy(cluster2) ;
2470 
2471  if(calorimeter==AliFiducialCut::kEMCAL)
2472  {
2473  GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
2474  GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
2475  }
2476 
2477  if(fPlotCluster)
2478  {
2479  //printf("Cells of cluster1: ");
2480  for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
2481  {
2482  //printf(" %d ",absIdList1[iDigit]);
2483 
2484  sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
2485 
2486  Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
2487  RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
2488 
2489  if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
2490  hCluster1->Fill(icol,irow,ecell*shareFraction1);
2491  else
2492  hCluster1->Fill(icol,irow,ecell);
2493  }
2494 
2495  //printf(" \n ");
2496  //printf("Cells of cluster2: ");
2497 
2498  for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
2499  {
2500  //printf(" %d ",absIdList2[iDigit]);
2501 
2502  sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
2503 
2504  Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
2505  RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
2506 
2507  if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) && absId2!=absIdList2[iDigit])
2508  hCluster2->Fill(icol,irow,ecell*shareFraction2);
2509  else
2510  hCluster2->Fill(icol,irow,ecell);
2511  }
2512  //printf(" \n ");
2513 
2514  gStyle->SetPadRightMargin(0.1);
2515  gStyle->SetPadLeftMargin(0.1);
2516  gStyle->SetOptStat(0);
2517  gStyle->SetOptFit(000000);
2518 
2519  if(maxCol-minCol > maxRow-minRow)
2520  {
2521  maxRow+= maxCol-minCol;
2522  }
2523  else
2524  {
2525  maxCol+= maxRow-minRow;
2526  }
2527 
2528  TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
2529  c->Divide(2,2);
2530  c->cd(1);
2531  gPad->SetGridy();
2532  gPad->SetGridx();
2533  gPad->SetLogz();
2534  hClusterMap ->SetAxisRange(minCol, maxCol,"X");
2535  hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
2536  hClusterMap ->Draw("colz TEXT");
2537  c->cd(2);
2538  gPad->SetGridy();
2539  gPad->SetGridx();
2540  gPad->SetLogz();
2541  hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
2542  hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
2543  hClusterLocMax ->Draw("colz TEXT");
2544  c->cd(3);
2545  gPad->SetGridy();
2546  gPad->SetGridx();
2547  gPad->SetLogz();
2548  hCluster1 ->SetAxisRange(minCol, maxCol,"X");
2549  hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
2550  hCluster1 ->Draw("colz TEXT");
2551  c->cd(4);
2552  gPad->SetGridy();
2553  gPad->SetGridx();
2554  gPad->SetLogz();
2555  hCluster2 ->SetAxisRange(minCol, maxCol,"X");
2556  hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
2557  hCluster2 ->Draw("colz TEXT");
2558 
2559  if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
2560  eventNumber,cluster->E(),nMax,ncells1,ncells2));
2561 
2562  delete c;
2563  delete hClusterMap;
2564  delete hClusterLocMax;
2565  delete hCluster1;
2566  delete hCluster2;
2567  }
2568 }
2569 
Bool_t fRecalculatePosition
Recalculate cluster position.
virtual Double_t Eta() const
AliEMCALGeometry * fEMCALGeo
! EMCAL geometry pointer.
Bool_t IsMCParticleInCalorimeterAcceptance(Int_t calo, AliVParticle *particle)
Check that a MC AOD is in the calorimeter acceptance.
double Double_t
Definition: External.C:58
void InitPHOSRecalibrationFactors()
Init PHOS recalibration factors.
Bool_t IsRecalculationOfClusterTrackMatchingOn() const
Definition: External.C:236
TGeoHMatrix * fEMCALMatrix[22]
Geometry matrices with alignments.
Int_t GetPHOSChannelStatus(Int_t imod, Int_t iCol, Int_t iRow) const
Float_t fLocMaxCutEDiff
Local maxima cut, when aggregating cells, next can be a bit higher.
const int debug
Definition: scanAll.C:15
Bool_t IsAbsIDsFromTCard(Int_t absId1, Int_t absId2, Int_t &rowDiff, Int_t &colDiff) const
TGeoHMatrix * fPHOSMatrix[5]
Geometry matrices with alignments.
Bool_t ClusterContainsBadChannel(Int_t calo, UShort_t *cellList, Int_t nCells)
AliPHOSGeoUtils * fPHOSGeo
! PHOS geometry pointer.
AliEMCALRecoUtils * GetEMCALRecoUtils() const
Float_t GetECross(Int_t absId, AliVCaloCells *cells, Int_t bc)
void SetEMCALL1PhaseInTimeRecalibrationForAllSM(TObjArray *map)
TString fileName
void SwitchOnDistToBadChannelRecalculation()
AliPHOSGeoUtils * GetPHOSGeometry() const
energy
Definition: HFPtSpectrum.C:44
Bool_t fOADBForEMCAL
Get calibration from OADB for EMCAL.
Float_t GetPHOSChannelRecalibrationFactor(Int_t imod, Int_t iCol, Int_t iRow) const
Float_t fLocMaxCutE
Local maxima cut must have more than this energy.
Float_t GetMCECellClusFracCorrection(Float_t eCell, Float_t eCluster) const
void SetEMCALChannelStatusMap(Int_t iSM, TH2I *h)
Bool_t fRunDependentCorrection
Switch on or off the recalibration dependent on T.
TCanvas * c
Definition: TestFitELoss.C:172
TObjArray * fPHOSRecalibrationFactors
Array of histograms with map of recalibration factors, PHOS.
TString fEMCALGeoName
Name of geometry to use for EMCAL.
TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const
Bool_t fRecalibration
Switch on or off the recalibration.
Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow) const
TH1F * GetEMCALChannelTimeRecalibrationFactors(Int_t bc) const
const TString calorimeter
Definition: anaM.C:36
Int_t fDebug
Debugging level.
TH1C * GetEMCALL1PhaseInTimeRecalibrationForAllSM() const
virtual UInt_t GetDetectorTag() const
Bool_t IsRecalibrationOn() const
virtual void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Int_t GetModuleNumber(AliCaloTrackParticle *particle, AliVEvent *inputEvent) const
Get the EMCAL/PHOS module number that corresponds to this particle.
Bool_t fLoadEMCALMatrices
Matrices set from configuration, not get from geometry.root or from ESDs/AODs.
AliEMCALGeometry * GetEMCALGeometry() const
int Int_t
Definition: External.C:63
Container for input particle information on CaloTrackCorr package.
Bool_t IsL1PhaseInTimeRecalibrationOn() const
Bool_t fRecalculateMatching
Recalculate cluster position.
Bool_t MaskFrameCluster(Int_t iSM, Int_t ieta) const
float Float_t
Definition: External.C:68
Bool_t fOADBForPHOS
Get calibration from OADB for PHOS.
Int_t * fMaskCellColumns
List of masked cells collumn index.
Bool_t fRemoveBadChannels
Check the channel status provided and remove clusters with bad channels.
void InitPHOSBadChannelStatusMap()
Init PHOS bad channels map.
void RecalculateClusterTrackMatching(AliVEvent *event, TObjArray *clusterArray=0x0, AliMCEvent *mc=0x0)
AliCalorimeterUtils()
Constructor. Initialize parameters.
virtual void InitParameters()
Initialize the parameters of the analysis.
Int_t fNMaskCellColumns
Number of masked columns.
Bool_t IsClusterSharedByTwoSuperModules(const AliEMCALGeometry *geom, AliVCluster *cluster)
Int_t GetNumberOfLocalMaxima(AliVCluster *cluster, AliVCaloCells *cells)
Find the number of local maxima in cluster.
Float_t fMCECellClusFracCorrParam[4]
Parameters for the function correcting the weight of the cells in the cluster.
Bool_t fEMCALGeoMatrixSet
Check if the transformation matrix is set for EMCAL.
TH2I * GetPHOSChannelStatusMap(Int_t imod) const
void SetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow, Double_t c=1)
TString GetPass()
Get passx from filename.
void RecalculateClusterPosition(AliVCaloCells *cells, AliVCluster *clu)
Bool_t fImportGeometryFromFile
Import geometry settings in geometry.root file.
Bool_t fCorrectELinearity
Correct cluster energy linearity.
Bool_t fPHOSGeoMatrixSet
Check if the transformation matrix is set for PHOS.
TString fOADBFilePathPHOS
Default path $ALICE_PHYSICS/OADB/PHOS, if needed change.
void CorrectClusterEnergy(AliVCluster *cl)
Correct cluster energy non linearity.
Float_t fCutZ
dZ cut on matching (EMCAL/PHOS).
Bool_t fLoadPHOSMatrices
Matrices set from configuration, not get from geometry.root or from ESDs/AODs.
void SetPHOSChannelStatusMap(Int_t imod, TH2I *h)
Int_t fNCellsFromPHOSBorder
Number of cells from PHOS border the cell with maximum amplitude has to be.
AliVTrack * GetMatchedTrack(AliVCluster *cluster, AliVEvent *event, Int_t index=-1) const
void RecalibrateCellTimeL1Phase(Double_t &time, Int_t calo, Int_t iSM, Int_t bunchCrossNumber) const
Recalculate time L1 phase shift if time recalibration available for EMCAL.
Bool_t fMCECellClusFracCorrOn
Correct or not the weight of cells in cluster.
Int_t fNSuperModulesUsed
Number of supermodules to be used in analysis, can be different than the real geo, to be used at initialization of histograms.
TString fImportGeometryFilePath
Path fo geometry.root file.
Bool_t fPlotCluster
Plot cluster in splitting method.
Int_t fRunNumber
Run number of the data, take it from data itself unless set by user.
virtual Double_t Phi() const
void RecalibrateCellAmplitude(Float_t &amp, Int_t calo, Int_t absId) const
Recalculate cell energy if recalibration factor.
void SetEMCALChannelTimeRecalibrationFactors(TObjArray *map)
Bool_t IsPHOSGeoMatrixSet() const
Bool_t IsTimeRecalibrationOn() const
Bool_t fOADBSet
AODB parameters already set.
unsigned short UShort_t
Definition: External.C:28
TObjArray * fPHOSBadChannelMap
Array of histograms with map of bad channels, PHOS.
TString fPHOSGeoName
Name of geometry to use for PHOS.
Bool_t AreNeighbours(Int_t calo, Int_t absId1, Int_t absId2) const
const char Option_t
Definition: External.C:48
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)
Int_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow) const
bool Bool_t
Definition: External.C:53
Bool_t IsEMCALGeoMatrixSet() const
Class with utils specific to calorimeter clusters/cells.
void AccessGeometry(AliVEvent *inputEvent)
AliEMCALRecoUtils * fEMCALRecoUtils
EMCAL utils for cluster rereconstruction.
void SetEMCALChannelRecalibrationFactors(Int_t iSM, TH2F *h)
void AccessOADB(AliVEvent *event)
TString fOADBFilePathEMCAL
Default path $ALICE_PHYSICS/OADB/EMCAL, if needed change.
void RecalibrateCellTime(Double_t &time, Int_t calo, Int_t absId, Int_t bunchCrossNumber) const
Recalculate time if time recalibration available for EMCAL not ready for PHOS.
Int_t GetMaxEnergyCell(AliVCaloCells *cells, AliVCluster *clu, Float_t &fraction) const
For a given CaloCluster, it gets the absId of the cell with maximum energy deposit.
Float_t fCutR
dR cut on matching (PHOS).
Float_t RecalibrateClusterEnergyWeightCell(AliVCluster *cluster, AliVCaloCells *cells, Float_t energyOrg)
virtual Int_t GetCaloLabel(Int_t i) const
void SplitEnergy(Int_t absId1, Int_t absId2, AliVCluster *cluster, AliVCaloCells *cells, AliAODCaloCluster *cluster1, AliAODCaloCluster *cluster2, Int_t nMax, Int_t eventNumber=0)
Bool_t GetFECCorrelatedCellAbsId(Int_t absId, Int_t absIdCorr[4]) const
Int_t GetModuleNumberCellIndexes(Int_t absId, Int_t calo, Int_t &icol, Int_t &irow, Int_t &iRCU) const
Get the EMCAL/PHOS module, columns, row and RCU/DDL number that corresponds to this absId...
TH2I * GetEMCALChannelStatusMap(Int_t iSM) const
void SetPHOSChannelRecalibrationFactor(Int_t imod, Int_t iCol, Int_t iRow, Double_t c=1)
Float_t RecalibrateClusterEnergy(AliVCluster *cluster, AliVCaloCells *cells)
Recalibrate the cluster energy, considering the recalibration map and the energy of the cells that co...
void GetEMCALSubregion(AliVCluster *clus, AliVCaloCells *cells, Int_t &regEta, Int_t &regPhi) const
Bool_t CheckCellFiducialRegion(AliVCluster *cluster, AliVCaloCells *cells) const
Int_t GetModuleNumberCellIndexesAbsCaloMap(Int_t absId, Int_t calo, Int_t &icol, Int_t &irow, Int_t &iRCU, Int_t &icolAbs, Int_t &irowAbs) const