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