AliPhysics  9c66e61 (9c66e61)
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"); // 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
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
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");
394  }
395  }
396  else
397  {
398  AliError("Do not calibrate L1 phase time");
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
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");
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
856  }
857 
858  // Recalibration factors
859 
860  if(bRecalE && ! bMC)
861  {
862  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on energy recalibration in EMCAL\n");
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");
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");
883  }
884 
885  // Recalculate position with method
886 
888 
889  // Non linearity
890 
891  if( bNonLin )
892  {
893  fCorrectELinearity = kTRUE;
894 
895  if(!bMC)
896  {
897  if ( debug > 0 )
898  printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kBeamTestCorrected xxx\n");
899 
901  }
902  else
903  {
904  if ( debug > 0 )
905  printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kPi0MCv3 xxx\n");
906 
908  }
909  } // Non linearity correction ON
910  else
911  {
912  fCorrectELinearity = kFALSE;
913 
914  if ( debug > 0 )
915  printf("ConfigureEMCALRecoUtils() xxx DON'T SET Non linearity correction xxx\n");
916 
918  } // Non linearity correction OFF
919 
920 }
921 
922 //_______________________________________________________________
924 //_______________________________________________________________
926 {
928 }
929 
930 //______________________________________________________________________________________
938 //______________________________________________________________________________________
939 Float_t AliCalorimeterUtils::GetECross(Int_t absID, AliVCaloCells* cells, Int_t bc)
940 {
941  if ( cells->IsEMCAL() )
942  {
943  Double_t tcell = cells->GetCellTime(absID);
944 
945  return fEMCALRecoUtils->GetECross(absID,tcell,cells,bc);
946  }
947  else // PHOS
948  {
949  Int_t icol = -1, irow = -1,iRCU = -1;
950  Int_t imod = GetModuleNumberCellIndexes(absID, AliFiducialCut::kPHOS, icol, irow, iRCU);
951 
952  Int_t absId1 = -1, absId2 = -1, absId3 = -1, absId4 = -1;
953 
954  Int_t relId1[] = { imod+1, 0, irow+1, icol };
955  Int_t relId2[] = { imod+1, 0, irow-1, icol };
956  Int_t relId3[] = { imod+1, 0, irow , icol+1 };
957  Int_t relId4[] = { imod+1, 0, irow , icol-1 };
958 
959  fPHOSGeo->RelToAbsNumbering(relId1, absId1);
960  fPHOSGeo->RelToAbsNumbering(relId2, absId2);
961  fPHOSGeo->RelToAbsNumbering(relId3, absId3);
962  fPHOSGeo->RelToAbsNumbering(relId4, absId4);
963 
964  Float_t ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0;
965 
966  if(absId1 > 0 ) ecell1 = cells->GetCellAmplitude(absId1);
967  if(absId2 > 0 ) ecell2 = cells->GetCellAmplitude(absId2);
968  if(absId3 > 0 ) ecell3 = cells->GetCellAmplitude(absId3);
969  if(absId4 > 0 ) ecell4 = cells->GetCellAmplitude(absId4);
970 
971  return ecell1+ecell2+ecell3+ecell4;
972  }
973 }
974 
975 //_______________________________________________________________
986 //______________________________________________________________________________
987 void AliCalorimeterUtils::GetEMCALSubregion(AliVCluster * clus, AliVCaloCells * cells,
988  Int_t & regEta, Int_t & regPhi) const
989 {
990  regEta = regPhi = -1 ;
991 
992  if(!clus->IsEMCAL()) return ;
993 
994  Int_t icol = -1, irow = -1, iRCU = -1;
995  Float_t clusterFraction = 0;
996 
997  Int_t absId = GetMaxEnergyCell(cells,clus,clusterFraction);
998 
999  Int_t sm = GetModuleNumberCellIndexes(absId,AliFiducialCut::kEMCAL,icol,irow,iRCU);
1000 
1001  // Shift by 48 to for impair SM
1002  if( sm%2 == 1) icol+=AliEMCALGeoParams::fgkEMCALCols;
1003 
1004  // Avoid borders
1005  if(icol < 2 || icol > 93 || irow < 2 || irow > 21) return;
1006 
1007  //
1008  // Eta regions
1009  //
1010 
1011  // Region 0: center of SM ~0.18<|eta|<0.55
1012  if ( icol > 9 && icol < 34 ) regEta = 0;
1013  else if ( icol > 62 && icol < 87 ) regEta = 0;
1014 
1015  // Region 3: frames ~0.1<|eta|<~0.22 ~0.51<|eta|<0.62
1016 
1017  else if ( icol <= 9 && icol >= 5 ) regEta = 3;
1018  else if ( icol <= 38 && icol >= 34 ) regEta = 3;
1019  else if ( icol <= 62 && icol >= 58 ) regEta = 3;
1020  else if ( icol <= 91 && icol >= 87 ) regEta = 3;
1021 
1022  // Region 1: |eta| < ~0.15
1023 
1024  else if ( icol < 58 && icol > 38 ) regEta = 1 ;
1025 
1026  // Region 2: |eta| > ~0.6
1027 
1028  else regEta = 2 ;
1029 
1030  //
1031  // Phi regions
1032  //
1033 
1034  if ( irow >= 2 && irow <= 5 ) regPhi = 0; // External
1035  else if ( irow >= 18 && irow <= 21 ) regPhi = 0; // External
1036  else if ( irow >= 6 && irow <= 9 ) regPhi = 1; // Mid
1037  else if ( irow >= 14 && irow <= 17 ) regPhi = 1; // Mid
1038  else regPhi = 2; //10-13 Central
1039 
1040 }
1041 
1042 //________________________________________________________________________________________
1049 //________________________________________________________________________________________
1051 {
1052  // Get SM number
1053  Int_t sm = fEMCALGeo->GetSuperModuleNumber(absId);
1054 
1055  // Get first absId of SM
1056  Int_t absIdSMMin = fEMCALGeo->GetAbsCellIdFromCellIndexes(sm,0,1); // for absIds, first is in col 1, module/tower ordering map ...
1057 
1058  // Get reference n and correlated 3
1059  for(Int_t k = 0; k < 4; k++ )
1060  {
1061  for(Int_t p = 0; p < 72; p++ )
1062  {
1063  Int_t n = absIdSMMin + 2*k + 16 *p;
1064 
1065  if ( absId == n || absId == n+1 ||
1066  absId == n+8 || absId == n+9 )
1067  {
1068  absIdCorr[0] = n ;
1069  absIdCorr[1] = n+1 ;
1070  absIdCorr[2] = n+8 ;
1071  absIdCorr[3] = n+9 ;
1072 
1073  //printf("n=%d, n+1=%d, n+8=%d, n+9=%d\n",
1074  // absIdCorr[0],absIdCorr[1],absIdCorr[2],absIdCorr[3]);
1075 
1076  return kTRUE;
1077  }
1078  }
1079  }
1080 
1081  // Not found;
1082  absIdCorr[0] = -1 ;
1083  absIdCorr[1] = -1 ;
1084  absIdCorr[2] = -1 ;
1085  absIdCorr[3] = -1 ;
1086 
1087  return kFALSE;
1088 }
1089 
1090 //________________________________________________________________________________________
1099 //________________________________________________________________________________________
1101  Int_t & rowDiff, Int_t & colDiff) const
1102 {
1103  rowDiff = -100;
1104  colDiff = -100;
1105 
1106  if(absId1 == absId2) return kFALSE;
1107 
1108  // Check if in same SM, if not for sure not same TCard
1109  Int_t sm1 = fEMCALGeo->GetSuperModuleNumber(absId1);
1110  Int_t sm2 = fEMCALGeo->GetSuperModuleNumber(absId2);
1111  if ( sm1 != sm2 ) return kFALSE ;
1112 
1113  // Get the column and row of each absId
1114  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1115 
1116  Int_t col1, row1;
1117  fEMCALGeo->GetCellIndex(absId1,sm1,iTower,iIphi,iIeta);
1118  fEMCALGeo->GetCellPhiEtaIndexInSModule(sm1,iTower,iIphi, iIeta,row1,col1);
1119 
1120  Int_t col2, row2;
1121  fEMCALGeo->GetCellIndex(absId2,sm2,iTower,iIphi,iIeta);
1122  fEMCALGeo->GetCellPhiEtaIndexInSModule(sm2,iTower,iIphi, iIeta,row2,col2);
1123 
1124  Int_t row0 = Int_t(row1-row1%8);
1125  Int_t col0 = Int_t(col1-col1%2);
1126 
1127  Int_t rowDiff0 = row2-row0;
1128  Int_t colDiff0 = col2-col0;
1129 
1130  rowDiff = row1-row2;
1131  colDiff = col1-col2;
1132 
1133  // TCard is made by 2x8 towers
1134  if ( colDiff0 >=0 && colDiff0 < 2 && rowDiff0 >=0 && rowDiff0 < 8 )
1135  {
1136 
1137 // printf("\t absId (%d,%d), sm %d; col (%d,%d), colDiff %d; row (%d,%d),rowDiff %d\n",
1138 // absId1 , absId2, sm1,
1139 // col1, col2, colDiff,
1140 // row1, row2, rowDiff);
1141  return kTRUE ;
1142  }
1143  else
1144  return kFALSE;
1145 }
1146 
1147 
1148 //________________________________________________________________________________________
1150 //________________________________________________________________________________________
1152  AliVCluster * clu,
1153  Float_t & clusterFraction) const
1154 {
1155  if( !clu || !cells )
1156  {
1157  AliInfo("Cluster or cells pointer is null!");
1158  return -1;
1159  }
1160 
1161  Double_t eMax =-1.;
1162  Double_t eTot = 0.;
1163  Double_t eCell =-1.;
1164  Float_t fraction = 1.;
1165  Float_t recalFactor = 1.;
1166  Int_t cellAbsId =-1 , absId =-1 ;
1167  Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
1168 
1170  if(clu->IsPHOS()) calo = AliFiducialCut::kPHOS ;
1171 
1172  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
1173  {
1174  cellAbsId = clu->GetCellAbsId(iDig);
1175 
1176  fraction = clu->GetCellAmplitudeFraction(iDig);
1177  if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
1178 
1179  iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
1180 
1181  if(IsRecalibrationOn())
1182  {
1183  if(calo == AliFiducialCut::kEMCAL)
1184  recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
1185  else
1186  recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
1187  }
1188 
1189  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
1190 
1191  if(eCell > eMax)
1192  {
1193  eMax = eCell;
1194  absId = cellAbsId;
1195  }
1196 
1197  eTot+=eCell;
1198 
1199  }// cell loop
1200 
1201  if(eTot > 0.1)
1202  clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
1203  else
1204  clusterFraction =1.;
1205 
1206  return absId;
1207 }
1208 
1209 //___________________________________________________________________________________
1213 //___________________________________________________________________________________
1214 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(AliVCluster* cluster,
1215  AliVEvent* event, Int_t index) const
1216 {
1217  AliVTrack *track = 0x0;
1218 
1219  //
1220  // EMCAL case only when matching is recalculated
1221  //
1222  if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
1223  {
1224  Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
1225  //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
1226 
1227  if(trackIndex < 0 )
1228  AliInfo(Form("Wrong track index %d, from recalculation", trackIndex));
1229  else
1230  track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
1231 
1232  return track ;
1233  }
1234 
1235  //
1236  // Normal case, get info from ESD or AOD
1237  //
1238 
1239  // No tracks matched
1240  if( cluster->GetNTracksMatched() < 1 ) return 0x0;
1241 
1242  // At least one match
1243  Int_t iTrack = 0; // only one match for AODs with index 0.
1244 
1245  // ESDs
1246  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
1247  {
1248  if( index >= 0 ) iTrack = index;
1249  else iTrack = ((AliESDCaloCluster*)cluster)->GetTracksMatched()->At(0); //cluster->GetTrackMatchedIndex();
1250 
1251  track = dynamic_cast<AliVTrack*> ( event->GetTrack(iTrack) );
1252  }
1253  else // AODs
1254  {
1255  if( index > 0 ) iTrack = index;
1256 
1257  track = dynamic_cast<AliVTrack*>( cluster->GetTrackMatched(iTrack) );
1258  }
1259 
1260  return track ;
1261 }
1262 
1268 {
1269  if( eCluster <= 0 || eCluster < eCell )
1270  {
1271  AliWarning(Form("Bad values eCell=%f, eCluster %f",eCell,eCluster));
1272  return 1;
1273  }
1274 
1275  Float_t frac = eCell / eCluster;
1276 
1277  Float_t correction = fMCECellClusFracCorrParam[0] +
1278  TMath::Exp( frac*fMCECellClusFracCorrParam[2]+fMCECellClusFracCorrParam[1] ) +
1279  fMCECellClusFracCorrParam[3]/TMath::Sqrt(frac);
1280 
1281 // printf("AliCalorimeterUtils::GetMCECellClusFracCorrection(eCell=%f, eCluster %f, frac %f) = %f\n",eCell, eCluster, frac, correction);
1282 // printf("\t %2.2f + TMath::Exp( %2.3f*%2.2f + %2.2f ) + %2.2f/TMath::Sqrt(%2.3f)) = %f\n",
1283 // fMCECellClusFracCorrParam[0],frac,fMCECellClusFracCorrParam[2],fMCECellClusFracCorrParam[1],fMCECellClusFracCorrParam[3], frac, correction);
1284 
1285  return correction;
1286 }
1287 
1288 //_____________________________________________________________________________________________________
1290 //_____________________________________________________________________________________________________
1291 Int_t AliCalorimeterUtils::GetModuleNumber(AliCaloTrackParticle * particle, AliVEvent * inputEvent) const
1292 {
1293  Int_t absId = -1;
1294 
1295  if(particle->GetDetectorTag()==AliFiducialCut::kEMCAL)
1296  {
1297  fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
1298 
1299  AliDebug(2,Form("EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d",
1300  particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId)));
1301 
1302  return fEMCALGeo->GetSuperModuleNumber(absId) ;
1303  } // EMCAL
1304  else if ( particle->GetDetectorTag() == AliFiducialCut::kPHOS )
1305  {
1306  // In case we use the MC reader, the input are TParticles,
1307  // in this case use the corresponing method in PHOS Geometry to get the particle.
1308  if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
1309  {
1310  Int_t mod =-1;
1311  Double_t z = 0., x=0.;
1312 
1313  TParticle* primary = ((AliMCEvent*)inputEvent)->Particle(particle->GetCaloLabel(0));
1314 
1315  if(primary)
1316  {
1317  fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
1318  }
1319  else
1320  {
1321  AliFatal("Primary not available, stop!");
1322  }
1323  return mod;
1324  }
1325  // Input are ESDs or AODs, get the PHOS module number like this.
1326  else
1327  {
1328  //FIXME
1329  //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
1330  //return GetModuleNumber(cluster);
1331  //MEFIX
1332  return -1;
1333  }
1334  } // PHOS
1335 
1336  return -1;
1337 }
1338 
1339 //_____________________________________________________________________
1341 //_____________________________________________________________________
1342 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
1343 {
1344  if(!cluster)
1345  {
1346  AliDebug(1,"AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
1347 
1348  return -1;
1349  }
1350 
1351  if ( cluster->GetNCells() <= 0 ) return -1;
1352 
1353  Int_t absId = cluster->GetCellAbsId(0);
1354 
1355  if ( absId < 0 ) return -1;
1356 
1357  if( cluster->IsEMCAL() )
1358  {
1359  AliDebug(2,Form("EMCAL absid %d, SuperModule %d",absId, fEMCALGeo->GetSuperModuleNumber(absId)));
1360 
1361  return fEMCALGeo->GetSuperModuleNumber(absId) ;
1362  } // EMCAL
1363  else if ( cluster->IsPHOS() )
1364  {
1365  Int_t relId[4];
1366  fPHOSGeo->AbsToRelNumbering(absId,relId);
1367 
1368  if (relId[1] != 0 ) return -1; // skip CPV only PHOS
1369 
1370  AliDebug(2,Form("PHOS absid %d Module %d",absId, relId[0]-1));
1371 
1372  return relId[0]-1;
1373  } // PHOS
1374 
1375  return -1;
1376 }
1377 
1378 //___________________________________________________________________________________________________
1380 //___________________________________________________________________________________________________
1382  Int_t & icol, Int_t & irow, Int_t & iRCU) const
1383 {
1384  Int_t imod = -1;
1385 
1386  if ( absId < 0 ) return -1 ;
1387 
1388  if ( calo == AliFiducialCut::kEMCAL || calo == AliFiducialCut::kDCAL)
1389  {
1390  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1391  fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1392  fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1393 
1394  if(imod < 0 || irow < 0 || icol < 0 )
1395  {
1396  AliFatal(Form("Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name",imod,icol,irow));
1397  }
1398 
1399  // In case of DCal C side, shift columns to match offline/online numbering
1400  // when calculating the DDL for Run2
1401  // See AliEMCALRawUtils::Digits2Raw and Raw2Digits.
1402  Int_t ico2 = icol ;
1403  if ( imod == 13 || imod == 15 || imod == 17 ) ico2 += 16;
1404 
1405  // RCU / DDL
1406  if(imod < 10 || (imod > 11 && imod < 18)) // (EMCAL Full || DCAL 2/3)
1407  {
1408  // RCU0 / DDL0
1409  if ( 0 <= irow && irow < 8 ) iRCU = 0; // first cable row
1410  else if ( 8 <= irow && irow < 16 &&
1411  0 <= ico2 && ico2 < 24 ) iRCU = 0; // first half;
1412  //second cable row
1413 
1414  // RCU1 / DDL1
1415  else if ( 8 <= irow && irow < 16 &&
1416  24 <= ico2 && ico2 < 48 ) iRCU = 1; // second half;
1417  //second cable row
1418  else if ( 16 <= irow && irow < 24 ) iRCU = 1; // third cable row
1419 
1420  if ( imod%2 == 1 ) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1421  }
1422  else
1423  {
1424  // 1/3 SM have one single SRU, just assign RCU/DDL 0
1425  iRCU = 0 ;
1426  }
1427 
1428  if ( iRCU < 0 )
1429  AliFatal(Form("Wrong EMCAL RCU number = %d", iRCU));
1430 
1431  return imod ;
1432  } // EMCAL
1433  else // PHOS
1434  {
1435  Int_t relId[4];
1436  fPHOSGeo->AbsToRelNumbering(absId,relId);
1437 
1438  if (relId[1] != 0 ) return -1; // skip CPV only PHOS
1439 
1440  irow = relId[2];
1441  icol = relId[3];
1442  imod = relId[0]-1;
1443  iRCU= (Int_t)(relId[2]-1)/16 ;
1444 
1445  //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
1446 
1447  if ( iRCU >= 4 )
1448  AliFatal(Form("Wrong PHOS RCU number = %d", iRCU));
1449 
1450  return imod;
1451  } // PHOS
1452 
1453  return -1;
1454 }
1455 
1456 //___________________________________________________________________________________________________
1459 //___________________________________________________________________________________________________
1461  Int_t & icol , Int_t & irow , Int_t & iRCU,
1462  Int_t & icolAbs, Int_t & irowAbs) const
1463 {
1464  Int_t imod = GetModuleNumberCellIndexes(absId, calo, icol, irow,iRCU);
1465 
1466  icolAbs = icol;
1467  irowAbs = irow;
1468 
1469  //
1470  // PHOS
1471  //
1472  if(calo == AliFiducialCut::kPHOS)
1473  {
1474  irowAbs = irow + 64 * imod;
1475 
1476  return imod;
1477  }
1478  //
1479  // EMCal/DCal
1480  //
1481  else
1482  {
1483  //
1484  // Shift collumns in even SM
1485  Int_t shiftEta = 48;
1486 
1487  // Shift collumn even more due to smaller acceptance of DCal collumns
1488  if ( imod > 11 && imod < 18) shiftEta+=48/3;
1489 
1490  icolAbs = (imod % 2) ? icol + shiftEta : icol;
1491 
1492  //
1493  // Shift rows per sector
1494  irowAbs = irow + 24 * Int_t(imod / 2);
1495 
1496  // Shift row less due to smaller acceptance of SM 10 and 11 to count DCal rows
1497  if ( imod > 11 && imod < 20) irowAbs -= (2*24 / 3);
1498 
1499  return imod ;
1500  }
1501 }
1502 
1503 
1504 //___________________________________________________________________________________________
1506 //___________________________________________________________________________________________
1507 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1508 {
1509  const Int_t nc = cluster->GetNCells();
1510 
1511  Int_t absIdList[nc];
1512  Float_t maxEList[nc];
1513 
1514  Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1515 
1516  return nMax;
1517 }
1518 
1519 //___________________________________________________________________________________________
1521 //___________________________________________________________________________________________
1522 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1523  Int_t *absIdList, Float_t *maxEList)
1524 {
1525  Int_t iDigitN = 0 ;
1526  Int_t iDigit = 0 ;
1527  Int_t absId1 = -1 ;
1528  Int_t absId2 = -1 ;
1529  const Int_t nCells = cluster->GetNCells();
1530 
1531  Float_t eCluster = RecalibrateClusterEnergy(cluster, cells);// recalculate cluster energy, avoid non lin correction.
1532 
1533  Float_t simuTotWeight = 0;
1535  {
1536  simuTotWeight = RecalibrateClusterEnergyWeightCell(cluster, cells,eCluster);// same but apply a weight
1537  simuTotWeight/= eCluster;
1538  }
1539 
1541  if(!cluster->IsEMCAL()) calorimeter = AliFiducialCut::kPHOS;
1542 
1543  //printf("cluster : ncells %d \n",nCells);
1544 
1545  Float_t emax = 0;
1546  Int_t idmax =-1;
1547  for(iDigit = 0; iDigit < nCells ; iDigit++)
1548  {
1549  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
1550  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1551  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1552 
1554  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1555 
1556  if( en > emax )
1557  {
1558  emax = en ;
1559  idmax = absIdList[iDigit] ;
1560  }
1561  //Int_t icol = -1, irow = -1, iRCU = -1;
1562  //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1563  //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
1564  }
1565 
1566  for(iDigit = 0 ; iDigit < nCells; iDigit++)
1567  {
1568  if( absIdList[iDigit] >= 0 )
1569  {
1570  absId1 = cluster->GetCellsAbsId()[iDigit];
1571 
1572  Float_t en1 = cells->GetCellAmplitude(absId1);
1573  RecalibrateCellAmplitude(en1,calorimeter,absId1);
1574 
1576  en1*=GetMCECellClusFracCorrection(en1,eCluster)/simuTotWeight;
1577 
1578  //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
1579 
1580  for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1581  {
1582  absId2 = cluster->GetCellsAbsId()[iDigitN] ;
1583 
1584  if(absId2==-1 || absId2==absId1) continue;
1585 
1586  //printf("\t %d : absIDj %d\n",iDigitN, absId2);
1587 
1588  Float_t en2 = cells->GetCellAmplitude(absId2);
1589  RecalibrateCellAmplitude(en2,calorimeter,absId2);
1590 
1592  en2*=GetMCECellClusFracCorrection(en2,eCluster)/simuTotWeight;
1593 
1594  //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
1595 
1596  if ( AreNeighbours(calorimeter, absId1, absId2) )
1597  {
1598  // printf("\t \t Neighbours \n");
1599  if ( en1 > en2 )
1600  {
1601  absIdList[iDigitN] = -1 ;
1602  //printf("\t \t indexN %d not local max\n",iDigitN);
1603  // but may be digit too is not local max ?
1604  if(en1 < en2 + fLocMaxCutEDiff) {
1605  //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
1606  absIdList[iDigit] = -1 ;
1607  }
1608  }
1609  else
1610  {
1611  absIdList[iDigit] = -1 ;
1612  //printf("\t \t index %d not local max\n",iDigitN);
1613  // but may be digitN too is not local max ?
1614  if(en1 > en2 - fLocMaxCutEDiff)
1615  {
1616  absIdList[iDigitN] = -1 ;
1617  //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
1618  }
1619  }
1620  } // if Are neighbours
1621  //else printf("\t \t NOT Neighbours \n");
1622  } // while digitN
1623  } // slot not empty
1624  } // while digit
1625 
1626  iDigitN = 0 ;
1627  for(iDigit = 0; iDigit < nCells; iDigit++)
1628  {
1629  if( absIdList[iDigit] >= 0 )
1630  {
1631  absIdList[iDigitN] = absIdList[iDigit] ;
1632 
1633  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1634  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1635 
1637  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1638 
1639  if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1640 
1641  maxEList[iDigitN] = en ;
1642 
1643  //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1644  iDigitN++ ;
1645  }
1646  }
1647 
1648  if ( iDigitN == 0 )
1649  {
1650  AliDebug(1,Form("No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f",
1651  idmax,emax,cluster->E()));
1652  iDigitN = 1 ;
1653  maxEList[0] = emax ;
1654  absIdList[0] = idmax ;
1655  }
1656 
1657 
1658  AliDebug(1,Form("In cluster E %2.2f (wth non lin. %2.2f), M02 %2.2f, M20 %2.2f, N maxima %d",
1659  cluster->E(),eCluster, cluster->GetM02(),cluster->GetM20(), iDigitN));
1660 
1661 // if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1662 // {
1663 // printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1664 // }
1665 
1666  return iDigitN ;
1667 }
1668 
1669 //____________________________________
1671 //____________________________________
1673 {
1674  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1675  {
1676  AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null");
1677  return TString("");
1678  }
1679 
1680  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1681  {
1682  AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null");
1683  return TString("");
1684  }
1685 
1686  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1687  if (pass.Contains("ass1")) return TString("pass1");
1688  else if (pass.Contains("ass2")) return TString("pass2");
1689  else if (pass.Contains("ass3")) return TString("pass3");
1690  else if (pass.Contains("ass4")) return TString("pass4");
1691  else if (pass.Contains("ass5")) return TString("pass5");
1692  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1693  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1694  {
1695  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1696  return TString("pass1");
1697  }
1698  else if (pass.Contains("LHC14a1a"))
1699  {
1700  AliInfo("Check that Energy calibration was enabled for this MC production in the tender, clusterizer or here!!");
1701 
1702  return TString("LHC14a1a");
1703  }
1704 
1705  // No condition fullfilled, give a default value
1706  AliInfo("Pass number string not found");
1707  return TString("");
1708 }
1709 
1710 //________________________________________
1712 //________________________________________
1714 {
1715  fEMCALGeoName = "";
1716  fPHOSGeoName = "";
1717 
1718  fEMCALGeoMatrixSet = kFALSE;
1719  fPHOSGeoMatrixSet = kFALSE;
1720 
1721  fRemoveBadChannels = kFALSE;
1722 
1724 
1725  fLocMaxCutE = 0.1 ;
1726  fLocMaxCutEDiff = 0.0 ;
1727 
1728  // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1729  // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1730  // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1731  // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1732  // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1733  // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
1734 
1735  fOADBSet = kFALSE;
1736  fOADBForEMCAL = kTRUE ;
1737  fOADBForPHOS = kFALSE;
1738 
1739  fOADBFilePathEMCAL = "" ; // $ALICE_PHYSICS/OADB/EMCAL
1740  fOADBFilePathPHOS = "$ALICE_PHYSICS/OADB/PHOS" ;
1741 
1742  fImportGeometryFromFile = kTRUE;
1744 
1745  fNSuperModulesUsed = 22;
1746 
1747  fMCECellClusFracCorrParam[0] = 0.78;
1748  fMCECellClusFracCorrParam[1] =-1.8;
1749  fMCECellClusFracCorrParam[2] =-6.3;
1750  fMCECellClusFracCorrParam[3] = 0.014;
1751 }
1752 
1753 //_____________________________________________________
1755 //_____________________________________________________
1757 {
1758  AliDebug(1,"Init bad channel map");
1759 
1760  // In order to avoid rewriting the same histograms
1761  Bool_t oldStatus = TH1::AddDirectoryStatus();
1762  TH1::AddDirectory(kFALSE);
1763 
1764  fPHOSBadChannelMap = new TObjArray(5);
1765  for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),
1766  Form("PHOS_BadMap_mod%d",i),
1767  64, 0, 64, 56, 0, 56));
1768 
1769  fPHOSBadChannelMap->SetOwner(kTRUE);
1770  fPHOSBadChannelMap->Compress();
1771 
1772  //In order to avoid rewriting the same histograms
1773  TH1::AddDirectory(oldStatus);
1774 }
1775 
1776 //______________________________________________________
1778 //______________________________________________________
1780 {
1781  AliDebug(1,"Init recalibration map");
1782 
1783  // In order to avoid rewriting the same histograms
1784  Bool_t oldStatus = TH1::AddDirectoryStatus();
1785  TH1::AddDirectory(kFALSE);
1786 
1788  for (int i = 0; i < 5; i++)
1789  {
1790  fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),
1791  Form("PHOSRecalFactors_Mod%d",i),
1792  64, 0, 64, 56, 0, 56));
1793  }
1794 
1795  // Init the histograms with 1
1796  for (Int_t m = 0; m < 5; m++)
1797  {
1798  for (Int_t i = 0; i < 56; i++)
1799  {
1800  for (Int_t j = 0; j < 64; j++)
1801  {
1803  }
1804  }
1805  }
1806 
1807  fPHOSRecalibrationFactors->SetOwner(kTRUE);
1808  fPHOSRecalibrationFactors->Compress();
1809 
1810  // In order to avoid rewriting the same histograms
1811  TH1::AddDirectory(oldStatus);
1812 }
1813 
1818 {
1819  if (fEMCALGeo) return;
1820 
1821  AliDebug(1,Form(" for run=%d",fRunNumber));
1822 
1823  if(fEMCALGeoName=="")
1824  {
1825  fEMCALGeo = AliEMCALGeometry::GetInstanceFromRunNumber(fRunNumber);
1826  AliInfo(Form("Get EMCAL geometry name to <%s> for run %d",fEMCALGeo->GetName(),fRunNumber));
1827  }
1828  else
1829  {
1830  fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
1831  AliInfo(Form("Set EMCAL geometry name to <%s>",fEMCALGeoName.Data()));
1832  }
1833 
1834  // Init geometry, I do not like much to do it like this ...
1835  if(fImportGeometryFromFile && !gGeoManager)
1836  {
1837  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1838  {
1839  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1840  if (fRunNumber < 140000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2010.root").data();
1841  else if(fRunNumber < 171000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2011.root").data();
1842  else if(fRunNumber < 198000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2012.root").data(); // 2012-2013
1843  else fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2015.root").data(); // >= 2015
1844  }
1845 
1846  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1847 
1848  TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1849  }
1850  else if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1851 }
1852 
1857 {
1858  if (fPHOSGeo) return;
1859 
1860  AliDebug(1,Form(" for run=%d",fRunNumber));
1861 
1862  if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1863 
1864  fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
1865 
1866  //if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1867 }
1868 
1869 //______________________________________________________________________________________________
1870 // Check that a MC ESD is in the calorimeter acceptance
1871 //______________________________________________________________________________________________
1873 {
1874  if(!particle || (calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS)) return kFALSE ;
1875 
1876  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
1878  {
1879  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1880  return kFALSE ;
1881  }
1882 
1883  if(calo == AliFiducialCut::kPHOS )
1884  {
1885  Int_t mod = 0 ;
1886  Double_t x = 0, z = 0 ;
1887  return GetPHOSGeometry()->ImpactOnEmc( particle, mod, z, x);
1888  }
1889  else if(calo == AliFiducialCut::kEMCAL)
1890  {
1891  Int_t absID = 0 ;
1892  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(),absID);
1893  if(ok)
1894  {
1895  Int_t icol = -1, irow = -1, iRCU = -1, status = 0;
1896  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1897  Bool_t bad = GetEMCALChannelStatus(nModule,icol,irow,status);
1898  if( bad ) ok = kFALSE;
1899  }
1900 
1901  return ok ;
1902  }
1903 
1904  return kFALSE ;
1905 }
1906 
1907 //______________________________________________________________________________________________________
1909 //______________________________________________________________________________________________________
1911 {
1912  if(!particle || (calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS)) return kFALSE ;
1913 
1914  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
1916  {
1917  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1918  return kFALSE ;
1919  }
1920 
1921  Float_t phi = particle->Phi();
1922  if(phi < 0) phi+=TMath::TwoPi();
1923 
1924  if(calo == AliFiducialCut::kPHOS )
1925  {
1926  Int_t mod = 0 ;
1927  Double_t x = 0, z = 0 ;
1928  Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1929  return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), phi, mod, z, x) ;
1930  }
1931  else if(calo == AliFiducialCut::kEMCAL)
1932  {
1933  Int_t absID = 0 ;
1934  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),phi,absID);
1935  if(ok)
1936  {
1937  Int_t icol = -1, irow = -1, iRCU = -1, status = 0;
1938  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1939  Bool_t bad = GetEMCALChannelStatus(nModule,icol,irow,status);
1940  if( bad ) ok = kFALSE;
1941  }
1942 
1943  return ok ;
1944  }
1945 
1946  return kFALSE ;
1947 }
1948 
1949 //______________________________________________________________________________________________________
1951 //______________________________________________________________________________________________________
1953 {
1954  if(!particle || (calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS)) return kFALSE ;
1955 
1956  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
1958  {
1959  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1960  return kFALSE ;
1961  }
1962 
1963  Float_t phi = particle->Phi();
1964  if(phi < 0) phi+=TMath::TwoPi();
1965 
1966  if(calo == AliFiducialCut::kPHOS )
1967  {
1968  Int_t mod = 0 ;
1969  Double_t x = 0, z = 0 ;
1970  Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1971  return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), phi, mod, z, x) ;
1972  }
1973  else if(calo == AliFiducialCut::kEMCAL)
1974  {
1975  Int_t absID = 0 ;
1976  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),phi,absID);
1977  if(ok)
1978  {
1979  Int_t icol = -1, irow = -1, iRCU = -1, status = 0;
1980  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1981  Bool_t bad = GetEMCALChannelStatus(nModule,icol,irow,status);
1982  if( bad ) ok = kFALSE;
1983  }
1984 
1985  return ok ;
1986  }
1987 
1988  return kFALSE ;
1989 }
1990 
1991 //_____________________________________________________________________________________________________
1992 // Check that a TLorentzVector is in the calorimeter acceptance, give the cell number where it hit.
1993 //_____________________________________________________________________________________________________
1995  Float_t phiOrg, Int_t & absID)
1996 {
1997  if(calo!=AliFiducialCut::kEMCAL && calo!=AliFiducialCut::kPHOS) return kFALSE ;
1998 
1999  if( (!IsPHOSGeoMatrixSet () && calo == AliFiducialCut::kPHOS ) ||
2001  {
2002  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
2003  return kFALSE ;
2004  }
2005 
2006  Float_t phi = phiOrg;
2007  if(phi < 0) phi+=TMath::TwoPi();
2008 
2009  if(calo == AliFiducialCut::kPHOS )
2010  {
2011  Int_t mod = 0 ;
2012  Double_t x = 0, z = 0 ;
2013  Double_t vtx[]={0,0,0} ;
2014  return GetPHOSGeometry()->ImpactOnEmc(vtx, theta, phi, mod, z, x) ;
2015  }
2016  else if(calo == AliFiducialCut::kEMCAL)
2017  {
2018  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta,phi,absID);
2019  if(ok)
2020  {
2021  Int_t icol = -1, irow = -1, iRCU = -1, status = 0;
2022  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
2023  Bool_t bad = GetEMCALChannelStatus(nModule,icol,irow,status);
2024  if( bad ) ok = kFALSE;
2025  }
2026 
2027  return ok ;
2028  }
2029 
2030  return kFALSE ;
2031 }
2032 
2033 //_______________________________________________________________________
2036 //_______________________________________________________________________
2038 {
2039  Int_t icol = ieta;
2040  if ( iSM%2 ) icol+=48; // Impair SM, shift index [0-47] to [48-96]
2041 
2043  {
2044  for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
2045  {
2046  if(icol==fMaskCellColumns[imask]) return kTRUE;
2047  }
2048  }
2049 
2050  return kFALSE;
2051 }
2052 
2053 //_________________________________________________________
2055 //_________________________________________________________
2057 {
2058  if(! opt)
2059  return;
2060 
2061  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
2062  printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
2063  printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
2065  if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
2066  printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
2067  printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
2068  printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
2069  printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
2070 
2071  printf("Recalibrate time? %d, With L1 phase run by run? %d\n",IsTimeRecalibrationOn(),IsL1PhaseInTimeRecalibrationOn());
2072 
2073  printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
2074  printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
2075 
2076  printf(" \n") ;
2077 }
2078 
2079 //_____________________________________________________________________________________________
2081 //_____________________________________________________________________________________________
2083 {
2084  Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
2085  Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
2086 
2087  if (IsRecalibrationOn())
2088  {
2089  if(calo == AliFiducialCut::kPHOS)
2090  {
2091  amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
2092  }
2093  else
2094  {
2095  amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
2096  }
2097  }
2098 }
2099 
2100 //____________________________________________________________________________________________________
2102 //____________________________________________________________________________________________________
2104 {
2106  {
2107  GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
2108  }
2109 }
2110 
2111 
2112 //____________________________________________________________________________________________________
2114 //____________________________________________________________________________________________________
2115 void AliCalorimeterUtils::RecalibrateCellTimeL1Phase(Double_t & time, Int_t calo, Int_t iSM, Int_t bunchCrossNumber) const
2116 {
2118  {
2119  GetEMCALRecoUtils()->RecalibrateCellTimeL1Phase(iSM, bunchCrossNumber, time);
2120  }
2121 }
2122 
2123 
2124 //__________________________________________________________________________
2126 //__________________________________________________________________________
2128  AliVCaloCells * cells)
2129 {
2130  if ( !IsRecalibrationOn() ) return cluster->E();
2131 
2132  // Initialize some used variables
2133  Float_t frac = 0., energy = 0.;
2134 
2135  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
2136  UShort_t * index = cluster->GetCellsAbsId() ;
2137  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
2138  Int_t ncells = cluster->GetNCells();
2139 
2141  if(cluster->IsPHOS()) calo = AliFiducialCut::kPHOS ;
2142 
2143  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
2144  for(Int_t icell = 0; icell < ncells; icell++)
2145  {
2146  Int_t absId = index[icell];
2147 
2148  frac = fraction[icell];
2149  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
2150 
2151  Float_t amp = cells->GetCellAmplitude(absId);
2152  RecalibrateCellAmplitude(amp,calo, absId);
2153 
2154  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy: before cal %f; after cal %f",
2155  calo,frac,cells->GetCellAmplitude(absId),amp));
2156 
2157  energy += amp*frac;
2158  }
2159 
2160  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
2161 
2162  return energy;
2163 }
2164 
2165 //_______________________________________________________________________________________________________
2168 //_______________________________________________________________________________________________________
2170  AliVCaloCells * cells, Float_t energyOrg)
2171 {
2172  // Initialize some used variables
2173  Float_t frac = 0., energy = 0.;
2174 
2175  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
2176  UShort_t * index = cluster->GetCellsAbsId() ;
2177  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
2178  Int_t ncells = cluster->GetNCells();
2179 
2181  if(cluster->IsPHOS()) calo = AliFiducialCut::kPHOS ;
2182 
2183  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
2184  for(Int_t icell = 0; icell < ncells; icell++)
2185  {
2186  Int_t absId = index[icell];
2187 
2188  frac = fraction[icell];
2189  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
2190 
2191  Float_t amp = cells->GetCellAmplitude(absId);
2192  if ( IsRecalibrationOn() ) RecalibrateCellAmplitude(amp,calo, absId);
2193 
2194  amp*=GetMCECellClusFracCorrection(amp,energyOrg);
2195 
2196  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy %f",
2197  calo,frac,cells->GetCellAmplitude(absId)));
2198 
2199  energy += amp*frac;
2200  }
2201 
2202  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
2203 
2204  return energy;
2205 }
2206 
2207 //__________________________________________________________________________________________
2210 //__________________________________________________________________________________________
2211 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
2212 {
2213  fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
2214 }
2215 
2216 //________________________________________________________________________________
2223 //________________________________________________________________________________
2225  TObjArray* clusterArray,
2226  AliMCEvent* mc)
2227 {
2228  if (!fRecalculateMatching) return ;
2229 
2230  fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo,mc) ;
2231 
2232  Float_t dZ = 2000;
2233  Float_t dR = 2000;
2234 
2235  Int_t nClusters = event->GetNumberOfCaloClusters();
2236  if(clusterArray) nClusters = clusterArray->GetEntriesFast();
2237 
2238  AliVCluster * clus = 0;
2239 
2240  for (Int_t iclus = 0; iclus < nClusters ; iclus++)
2241  {
2242  if ( clusterArray ) clus = (AliVCluster*) clusterArray->At(iclus) ;
2243  else clus = event->GetCaloCluster(iclus) ;
2244 
2245  if (!clus->IsEMCAL()) continue ;
2246 
2247  //
2248  // Put track residuals in cluster
2249  //
2250  fEMCALRecoUtils->GetMatchedResiduals(clus->GetID(),dZ,dR);
2251 
2252  if ( TMath::Abs(clus->GetTrackDx()) < 500 )
2253  AliDebug(2,Form("Residuals (Old, New): z (%2.4f,%2.4f), x (%2.4f,%2.4f)",
2254  clus->GetTrackDz(),dZ,clus->GetTrackDx(),dR));
2255 
2256  clus->SetTrackDistance(dR,dZ);
2257 
2258  //
2259  // Remove old matches in cluster
2260  //
2261  if(clus->GetNTracksMatched() > 0)
2262  {
2263  if(!strcmp("AliESDCaloCluster",Form("%s",clus->ClassName())))
2264  {
2265  TArrayI arrayTrackMatched(0);
2266  ((AliESDCaloCluster*)clus)->AddTracksMatched(arrayTrackMatched);
2267  }
2268  else
2269  {
2270  for(Int_t iTrack = 0; iTrack < clus->GetNTracksMatched(); iTrack++)
2271  {
2272  ((AliAODCaloCluster*)clus)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)clus)->GetTrackMatched(iTrack));
2273  }
2274  }
2275  }
2276 
2277  //
2278  // Now put first track index in cluster.
2279  //
2280  Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(iclus);
2281  if ( trackIndex >= 0 )
2282  {
2283  if(!strcmp("AliESDCaloCluster",Form("%s",clus->ClassName())))
2284  {
2285  TArrayI arrayTrackMatched(1);
2286  arrayTrackMatched[0] = trackIndex;
2287  ((AliESDCaloCluster*)clus)->AddTracksMatched(arrayTrackMatched);
2288  }
2289  else
2290  {
2291  ((AliAODCaloCluster*)clus)->AddTrackMatched((TObject*)event->GetTrack(trackIndex));
2292  }
2293  }
2294 
2295  } // cluster loop
2296 }
2297 
2298 //___________________________________________________________________________
2311 //___________________________________________________________________________
2313  AliVCluster* cluster,
2314  AliVCaloCells* cells,
2315  //Float_t & e1, Float_t & e2,
2316  AliAODCaloCluster* cluster1,
2317  AliAODCaloCluster* cluster2,
2318  Int_t nMax, Int_t eventNumber)
2319 {
2320  TH2F* hClusterMap = 0 ;
2321  TH2F* hClusterLocMax = 0 ;
2322  TH2F* hCluster1 = 0 ;
2323  TH2F* hCluster2 = 0 ;
2324 
2325  if(fPlotCluster)
2326  {
2327  hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
2328  hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
2329  hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
2330  hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
2331  hClusterMap ->SetXTitle("column");
2332  hClusterMap ->SetYTitle("row");
2333  hClusterLocMax ->SetXTitle("column");
2334  hClusterLocMax ->SetYTitle("row");
2335  hCluster1 ->SetXTitle("column");
2336  hCluster1 ->SetYTitle("row");
2337  hCluster2 ->SetXTitle("column");
2338  hCluster2 ->SetYTitle("row");
2339  }
2340 
2342  if(cluster->IsPHOS())
2343  {
2344  calorimeter = AliFiducialCut::kPHOS;
2345  AliWarning("Not supported for PHOS yet");
2346  return;
2347  }
2348 
2349  const Int_t ncells = cluster->GetNCells();
2350  Int_t absIdList[ncells];
2351 
2352  Float_t e1 = 0, e2 = 0 ;
2353  Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
2354  Float_t eCluster = 0;
2355  Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
2356  for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
2357  {
2358  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
2359 
2360  Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
2361  RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
2362  eCluster+=ec;
2363 
2364  if(fPlotCluster)
2365  {
2366  //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
2367  sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
2368  if(sm > -1 && sm < 12) // just to avoid compilation warning
2369  {
2370  if(icol > maxCol) maxCol = icol;
2371  if(icol < minCol) minCol = icol;
2372  if(irow > maxRow) maxRow = irow;
2373  if(irow < minRow) minRow = irow;
2374  hClusterMap->Fill(icol,irow,ec);
2375  }
2376  }
2377  }
2378 
2379  // Init counters and variables
2380  Int_t ncells1 = 1 ;
2381  UShort_t absIdList1[9] ;
2382  Double_t fracList1 [9] ;
2383  absIdList1[0] = absId1 ;
2384  fracList1 [0] = 1. ;
2385 
2386  Float_t ecell1 = cells->GetCellAmplitude(absId1);
2387  RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
2388  e1 = ecell1;
2389 
2390  Int_t ncells2 = 1 ;
2391  UShort_t absIdList2[9] ;
2392  Double_t fracList2 [9] ;
2393  absIdList2[0] = absId2 ;
2394  fracList2 [0] = 1. ;
2395 
2396  Float_t ecell2 = cells->GetCellAmplitude(absId2);
2397  RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
2398  e2 = ecell2;
2399 
2400  if(fPlotCluster)
2401  {
2402  Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
2403  sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
2404  hClusterLocMax->Fill(icol1,irow1,ecell1);
2405  sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
2406  hClusterLocMax->Fill(icol2,irow2,ecell2);
2407  }
2408 
2409  // Very rough way to share the cluster energy
2410  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
2411  Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
2412  Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
2413 
2414  for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
2415  {
2416  Int_t absId = absIdList[iDigit];
2417 
2418  if ( absId==absId1 || absId==absId2 || absId < 0 ) continue;
2419 
2420  Float_t ecell = cells->GetCellAmplitude(absId);
2421  RecalibrateCellAmplitude(ecell, calorimeter, absId);
2422 
2423  if(AreNeighbours(calorimeter, absId1,absId ))
2424  {
2425  absIdList1[ncells1]= absId;
2426 
2427  if(AreNeighbours(calorimeter, absId2,absId ))
2428  {
2429  fracList1[ncells1] = shareFraction1;
2430  e1 += ecell*shareFraction1;
2431  }
2432  else
2433  {
2434  fracList1[ncells1] = 1.;
2435  e1 += ecell;
2436  }
2437 
2438  ncells1++;
2439 
2440  } // neigbour to cell1
2441 
2442  if(AreNeighbours(calorimeter, absId2,absId ))
2443  {
2444  absIdList2[ncells2]= absId;
2445 
2446  if(AreNeighbours(calorimeter, absId1,absId ))
2447  {
2448  fracList2[ncells2] = shareFraction2;
2449  e2 += ecell*shareFraction2;
2450  }
2451  else
2452  {
2453  fracList2[ncells2] = 1.;
2454  e2 += ecell;
2455  }
2456 
2457  ncells2++;
2458 
2459  } // neigbour to cell2
2460  }
2461 
2462  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",
2463  nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2));
2464 
2465  cluster1->SetE(e1);
2466  cluster2->SetE(e2);
2467 
2468  cluster1->SetNCells(ncells1);
2469  cluster2->SetNCells(ncells2);
2470 
2471  cluster1->SetCellsAbsId(absIdList1);
2472  cluster2->SetCellsAbsId(absIdList2);
2473 
2474  cluster1->SetCellsAmplitudeFraction(fracList1);
2475  cluster2->SetCellsAmplitudeFraction(fracList2);
2476 
2477  // Correct linearity
2478  CorrectClusterEnergy(cluster1) ;
2479  CorrectClusterEnergy(cluster2) ;
2480 
2481  if(calorimeter==AliFiducialCut::kEMCAL)
2482  {
2485  }
2486 
2487  if(fPlotCluster)
2488  {
2489  //printf("Cells of cluster1: ");
2490  for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
2491  {
2492  //printf(" %d ",absIdList1[iDigit]);
2493 
2494  sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
2495 
2496  Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
2497  RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
2498 
2499  if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
2500  hCluster1->Fill(icol,irow,ecell*shareFraction1);
2501  else
2502  hCluster1->Fill(icol,irow,ecell);
2503  }
2504 
2505  //printf(" \n ");
2506  //printf("Cells of cluster2: ");
2507 
2508  for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
2509  {
2510  //printf(" %d ",absIdList2[iDigit]);
2511 
2512  sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
2513 
2514  Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
2515  RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
2516 
2517  if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) && absId2!=absIdList2[iDigit])
2518  hCluster2->Fill(icol,irow,ecell*shareFraction2);
2519  else
2520  hCluster2->Fill(icol,irow,ecell);
2521  }
2522  //printf(" \n ");
2523 
2524  gStyle->SetPadRightMargin(0.1);
2525  gStyle->SetPadLeftMargin(0.1);
2526  gStyle->SetOptStat(0);
2527  gStyle->SetOptFit(000000);
2528 
2529  if(maxCol-minCol > maxRow-minRow)
2530  {
2531  maxRow+= maxCol-minCol;
2532  }
2533  else
2534  {
2535  maxCol+= maxRow-minRow;
2536  }
2537 
2538  TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
2539  c->Divide(2,2);
2540  c->cd(1);
2541  gPad->SetGridy();
2542  gPad->SetGridx();
2543  gPad->SetLogz();
2544  hClusterMap ->SetAxisRange(minCol, maxCol,"X");
2545  hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
2546  hClusterMap ->Draw("colz TEXT");
2547  c->cd(2);
2548  gPad->SetGridy();
2549  gPad->SetGridx();
2550  gPad->SetLogz();
2551  hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
2552  hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
2553  hClusterLocMax ->Draw("colz TEXT");
2554  c->cd(3);
2555  gPad->SetGridy();
2556  gPad->SetGridx();
2557  gPad->SetLogz();
2558  hCluster1 ->SetAxisRange(minCol, maxCol,"X");
2559  hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
2560  hCluster1 ->Draw("colz TEXT");
2561  c->cd(4);
2562  gPad->SetGridy();
2563  gPad->SetGridx();
2564  gPad->SetLogz();
2565  hCluster2 ->SetAxisRange(minCol, maxCol,"X");
2566  hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
2567  hCluster2 ->Draw("colz TEXT");
2568 
2569  if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
2570  eventNumber,cluster->E(),nMax,ncells1,ncells2));
2571 
2572  delete c;
2573  delete hClusterMap;
2574  delete hClusterLocMax;
2575  delete hCluster1;
2576  delete hCluster2;
2577  }
2578 }
2579 
Bool_t fRecalculatePosition
Recalculate cluster position.
virtual Double_t Eta() const
AliEMCALGeometry * fEMCALGeo
! EMCAL geometry pointer.
Bool_t IsL1PhaseInTimeRecalibrationOn() const
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.
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)
Int_t GetNumberOfCellsFromEMCALBorder() const
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
void SetNonLinearityFunction(Int_t fun)
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
Some utilities for cluster and cell treatment.
virtual void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
TH2I * GetEMCALChannelStatusMap(Int_t iSM) const
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.
Int_t GetMatchedTrackIndex(Int_t clsIndex)
AliEMCALGeometry * GetEMCALGeometry() const
void SwitchOnL1PhaseInTimeRecalibration()
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 IsEMCALNoBorderAtEta0() const
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)
Bool_t IsTimeRecalibrationOn() const
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)
void GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
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.
void SetPositionAlgorithm(Int_t alg)
void SwitchOffL1PhaseInTimeRecalibration()
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)
void SwitchOnDistToBadChannelRecalculation()
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).
void RecalculateClusterPosition(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
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.
void SwitchOnRejectExoticCluster()
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 SetExoticCellMinAmplitudeCut(Float_t ma)
void RecalibrateCellTimeL1Phase(Int_t iSM, Int_t bc, Double_t &time) const
void SetEMCALChannelTimeRecalibrationFactors(TObjArray *map)
Bool_t IsPHOSGeoMatrixSet() const
void FindMatches(AliVEvent *event, TObjArray *clusterArr=0x0, const AliEMCALGeometry *geom=0x0, AliMCEvent *mc=0x0)
Bool_t IsTimeRecalibrationOn() const
Bool_t fOADBSet
AODB parameters already set.
unsigned short UShort_t
Definition: External.C:28
void SetExoticCellFractionCut(Float_t f)
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)
bool Bool_t
Definition: External.C:53
Bool_t IsEMCALGeoMatrixSet() const
Class with utils specific to calorimeter clusters/cells.
Float_t CorrectClusterEnergyLinearity(AliVCluster *clu)
void AccessGeometry(AliVEvent *inputEvent)
AliEMCALRecoUtils * fEMCALRecoUtils
EMCAL utils for cluster rereconstruction.
Bool_t ClusterContainsBadChannel(const AliEMCALGeometry *geom, const UShort_t *cellList, Int_t nCells)
void SetEMCALChannelRecalibrationFactors(Int_t iSM, TH2F *h)
void AccessOADB(AliVEvent *event)
Float_t GetECross(Int_t absID, Double_t tcell, AliVCaloCells *cells, Int_t bc)
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).
Int_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow, Int_t status) const
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...
void RecalibrateCellTime(Int_t absId, Int_t bc, Double_t &time, Bool_t isLGon=kFALSE) const
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