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