AliPhysics  a0db429 (a0db429)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
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 
25 // --- ANALYSIS system ---
26 #include "AliCalorimeterUtils.h"
27 #include "AliESDEvent.h"
28 #include "AliMCEvent.h"
29 #include "AliStack.h"
30 #include "AliAODPWG4Particle.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 "AliLog.h"
38 
39 // --- Detector ---
40 #include "AliEMCALGeometry.h"
41 #include "AliPHOSGeoUtils.h"
42 
46 
47 
48 //____________________________________________
50 //____________________________________________
52 TObject(), fDebug(0),
53 fEMCALGeoName(""),
54 fPHOSGeoName (""),
55 fEMCALGeo(0x0), fPHOSGeo(0x0),
56 fEMCALGeoMatrixSet(kFALSE), fPHOSGeoMatrixSet(kFALSE),
57 fLoadEMCALMatrices(kFALSE), fLoadPHOSMatrices(kFALSE),
58 fRemoveBadChannels(kFALSE), fPHOSBadChannelMap(0x0),
59 fNCellsFromPHOSBorder(0),
60 fNMaskCellColumns(0), fMaskCellColumns(0x0),
61 fRecalibration(kFALSE), fRunDependentCorrection(kFALSE),
62 fPHOSRecalibrationFactors(), fEMCALRecoUtils(new AliEMCALRecoUtils),
63 fRecalculatePosition(kFALSE), fCorrectELinearity(kFALSE),
64 fRecalculateMatching(kFALSE),
65 fCutR(20), fCutZ(20),
66 fCutEta(20), fCutPhi(20),
67 fLocMaxCutE(0), fLocMaxCutEDiff(0),
68 fPlotCluster(0), fOADBSet(kFALSE),
69 fOADBForEMCAL(kFALSE), fOADBForPHOS(kFALSE),
70 fOADBFilePathEMCAL(""), fOADBFilePathPHOS(""),
71 fImportGeometryFromFile(0), fImportGeometryFilePath(""),
72 fNSuperModulesUsed(0), fRunNumber(0),
73 fMCECellClusFracCorrOn(0), fMCECellClusFracCorrParam()
74 {
76  for(Int_t i = 0; i < 22; i++) fEMCALMatrix[i] = 0 ;
77  for(Int_t i = 0; i < 5 ; i++) fPHOSMatrix [i] = 0 ;
78 }
79 
80 //_________________________________________
81 // Destructor.
82 //_________________________________________
84 {
85  //if(fPHOSGeo) delete fPHOSGeo ;
86  if(fEMCALGeo) delete fEMCALGeo ;
87 
89 {
90  fPHOSBadChannelMap->Clear();
91  delete fPHOSBadChannelMap;
92  }
93 
95 {
98  }
99 
100  if(fEMCALRecoUtils) delete fEMCALRecoUtils ;
101  if(fNMaskCellColumns) delete [] fMaskCellColumns;
102 }
103 
104 //____________________________________________________
107 //____________________________________________________
108 void AliCalorimeterUtils::AccessOADB(AliVEvent* event)
109 {
110  // Set it only once
111  if(fOADBSet) return ;
112 
113  if(fRunNumber <= 0) fRunNumber = event->GetRunNumber() ; // take the run number from the event itself
114  TString pass = GetPass();
115 
116  // EMCAL
117  if(fOADBForEMCAL)
118  {
119  AliInfo(Form("Get AODB parameters from EMCAL in %s for run %d, and <%s>",fOADBFilePathEMCAL.Data(),fRunNumber,pass.Data()));
120 
121  Int_t nSM = fEMCALGeo->GetNumberOfSuperModules();
122 
123  // Bad map
125  {
126  AliOADBContainer *contBC=new AliOADBContainer("");
127  contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePathEMCAL.Data()),"AliEMCALBadChannels");
128 
129  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRunNumber);
130 
131  if(arrayBC)
132  {
134  AliInfo("Remove EMCAL bad cells");
135 
136  for (Int_t i=0; i<nSM; ++i)
137  {
138  TH2I *hbm = GetEMCALChannelStatusMap(i);
139 
140  if (hbm)
141  delete hbm;
142 
143  hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
144 
145  if (!hbm)
146  {
147  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
148  continue;
149  }
150 
151  hbm->SetDirectory(0);
153 
154  } // loop
155  } else AliInfo("Do NOT remove EMCAL bad channels\n"); // run array
156 
157  delete contBC;
158  } // Remove bad
159 
160  // Energy Recalibration
161  if(fRecalibration)
162  {
163  AliOADBContainer *contRF=new AliOADBContainer("");
164 
165  contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRecalib");
166 
167  TObjArray *recal=(TObjArray*)contRF->GetObject(fRunNumber);
168 
169  if(recal)
170  {
171  TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
172 
173  if(recalpass)
174  {
175  TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
176 
177  if(recalib)
178  {
179  AliInfo("Recalibrate EMCAL");
180  for (Int_t i=0; i < nSM; ++i)
181  {
183 
184  if (h)
185  delete h;
186 
187  h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
188 
189  if (!h)
190  {
191  AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
192  continue;
193  }
194 
195  h->SetDirectory(0);
196 
198  } // SM loop
199  } else AliInfo("Do NOT recalibrate EMCAL, no params object array"); // array ok
200  } else AliInfo("Do NOT recalibrate EMCAL, no params for pass"); // array pass ok
201  } else AliInfo("Do NOT recalibrate EMCAL, no params for run"); // run number array ok
202 
203  delete contRF;
204  // once set, apply run dependent corrections if requested
205  //fEMCALRecoUtils->SetRunDependentCorrections(fRunNumber);
206 
207  } // Recalibration on
208 
209  // Energy Recalibration, apply on top of previous calibration factors
211  {
212  AliOADBContainer *contRFTD=new AliOADBContainer("");
213 
214  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALRunDepTempCalibCorrections");
215 
216  TH1S *htd=(TH1S*)contRFTD->GetObject(fRunNumber);
217 
218  //If it did not exist for this run, get closes one
219  if (!htd)
220  {
221  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",fRunNumber));
222  // let's get the closest fRunNumber instead then..
223  Int_t lower = 0;
224  Int_t ic = 0;
225  Int_t maxEntry = contRFTD->GetNumberOfEntries();
226 
227  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < fRunNumber) )
228  {
229  lower = ic;
230  ic++;
231  }
232 
233  Int_t closest = lower;
234  if ( (ic<maxEntry) &&
235  (contRFTD->LowerLimit(ic)-fRunNumber) < (fRunNumber - contRFTD->UpperLimit(lower)) )
236  {
237  closest = ic;
238  }
239 
240  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
241  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
242  }
243 
244  if(htd)
245  {
246  AliInfo("Recalibrate (Temperature) EMCAL");
247 
248  for (Int_t ism=0; ism<nSM; ++ism)
249  {
250  for (Int_t icol=0; icol<48; ++icol)
251  {
252  for (Int_t irow=0; irow<24; ++irow)
253  {
254  Float_t factor = GetEMCALChannelRecalibrationFactor(ism,icol,irow);
255 
256  Int_t absID = fEMCALGeo->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
257  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
258 
259  //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
260  // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
261 
262  SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
263  } // columns
264  } // rows
265  } // SM loop
266  } else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
267 
268  delete contRFTD;
269  } // Run by Run T calibration
270 
271  // Time Recalibration
272  if(fEMCALRecoUtils->IsTimeRecalibrationOn())
273  {
274  AliOADBContainer *contTRF=new AliOADBContainer("");
275 
276  contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePathEMCAL.Data()),"AliEMCALTimeCalib");
277 
278  TObjArray *trecal=(TObjArray*)contTRF->GetObject(fRunNumber);
279 
280  if(trecal)
281  {
282  TString passM = pass;
283  if(pass=="spc_calo") passM = "pass3";
284  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
285 
286  if(trecalpass)
287  {
288  AliInfo("Time Recalibrate EMCAL");
289  for (Int_t ibc = 0; ibc < 4; ++ibc)
290  {
292 
293  if (h)
294  delete h;
295 
296  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
297 
298  if (!h)
299  {
300  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
301  continue;
302  }
303 
304  h->SetDirectory(0);
305 
307  } // bunch crossing loop
308  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
309  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
310 
311  delete contTRF;
312  } // Recalibration on
313 
314  }// EMCAL
315 
316  // PHOS
317  if(fOADBForPHOS)
318  {
319  AliInfo(Form("Get AODB parameters from PHOS in %s for run %d, and <%s>",fOADBFilePathPHOS.Data(),fRunNumber,pass.Data()));
320 
321  // Bad map
323  {
324  AliOADBContainer badmapContainer(Form("phosBadMap"));
325  TString fileName="$ALICE_PHYSICS/OADB/PHOS/PHOSBadMaps.root";
326  badmapContainer.InitFromFile(Form("%s/PHOSBadMaps.root",fOADBFilePathPHOS.Data()),"phosBadMap");
327 
328  //Use a fixed run number from year 2010, this year not available yet.
329  TObjArray *maps = (TObjArray*)badmapContainer.GetObject(139000,"phosBadMap");
330  if(!maps)
331  {
332  AliInfo(Form("Can not read PHOS bad map for run %d",fRunNumber)) ;
333  }
334  else
335  {
336  AliInfo(Form("Setting PHOS bad map with name %s",maps->GetName())) ;
337 
338  for(Int_t mod = 1; mod < 5; mod++)
339  {
340  TH2I *hbmPH = GetPHOSChannelStatusMap(mod);
341 
342  if(hbmPH)
343  delete hbmPH ;
344 
345  hbmPH = (TH2I*)maps->At(mod);
346 
347  if(hbmPH) hbmPH->SetDirectory(0);
348 
349  SetPHOSChannelStatusMap(mod-1,hbmPH);
350 
351  } // modules loop
352  } // maps exist
353  } // Remove bad channels
354  } // PHOS
355 
356  // Parameters already set once, so do not it again
357  fOADBSet = kTRUE;
358 }
359 
360 //_____________________________________________________________
363 //_____________________________________________________________
364 void AliCalorimeterUtils::AccessGeometry(AliVEvent* inputEvent)
365 {
366  // First init the geometry, a priory not done before
367  if(fRunNumber <=0 ) fRunNumber = inputEvent->GetRunNumber() ;
368 
369  InitPHOSGeometry ();
371 
372  //Get the EMCAL transformation geometry matrices from ESD
374  {
376  {
377  AliInfo("Load user defined EMCAL geometry matrices");
378 
379  // OADB if available
380  AliOADBContainer emcGeoMat("AliEMCALgeo");
381  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePathEMCAL.Data()),"AliEMCALgeo");
382  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(fRunNumber,"EmcalMatrices");
383 
384  for(Int_t mod=0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
385  {
386  if (!fEMCALMatrix[mod]) // Get it from OADB
387  {
388  AliDebug(1,Form("EMCAL matrices SM %d, %p", mod,((TGeoHMatrix*) matEMCAL->At(mod))));
389  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
390 
391  fEMCALMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
392  }
393 
394  if(fEMCALMatrix[mod])
395  {
396  if(fDebug > 1)
397  fEMCALMatrix[mod]->Print();
398 
399  fEMCALGeo->SetMisalMatrix(fEMCALMatrix[mod],mod) ;
400  }
401  else if(gGeoManager)
402  {
403  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
404  fEMCALGeo->SetMisalMatrix(fEMCALGeo->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
405  }
406  else
407  {
408  AliError(Form("Alignment atrix for SM %d is not available",mod));
409  }
410  } // SM loop
411 
412  fEMCALGeoMatrixSet = kTRUE;//At least one, so good
413 
414  } // Load matrices
415  else if (!gGeoManager)
416  {
417  AliDebug(1,"Load EMCAL misalignment matrices");
418  if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
419  {
420  for(Int_t mod = 0; mod < (fEMCALGeo->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
421  {
422  //printf("Load ESD matrix %d, %p\n",mod,((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod));
423  if(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod))
424  {
425  fEMCALGeo->SetMisalMatrix(((AliESDEvent*)inputEvent)->GetEMCALMatrix(mod),mod) ;
426  }
427  }// loop over super modules
428 
429  fEMCALGeoMatrixSet = kTRUE;//At least one, so good
430 
431  }//ESD as input
432  else
433  {
434  AliDebug(1,"Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file");
435  }//AOD as input
436  }//Get matrix from data
437  else if(gGeoManager)
438  {
439  fEMCALGeoMatrixSet = kTRUE;
440  }
441  }//EMCAL geo && no geoManager
442 
443  //Get the PHOS transformation geometry matrices from ESD
445  {
447  {
448  AliInfo("Load user defined PHOS geometry matrices");
449 
450  // OADB if available
451  AliOADBContainer geomContainer("phosGeo");
452  geomContainer.InitFromFile(Form("%s/PHOSGeometry.root",fOADBFilePathPHOS.Data()),"PHOSRotationMatrixes");
453  TObjArray *matPHOS = (TObjArray*)geomContainer.GetObject(139000,"PHOSRotationMatrixes");
454 
455  for(Int_t mod = 0 ; mod < 5 ; mod++)
456  {
457  if (!fPHOSMatrix[mod]) // Get it from OADB
458  {
459  AliDebug(1,Form("PHOS matrices module %d, %p",mod,((TGeoHMatrix*) matPHOS->At(mod))));
460  //((TGeoHMatrix*) matPHOS->At(mod))->Print();
461 
462  fPHOSMatrix[mod] = (TGeoHMatrix*) matPHOS->At(mod) ;
463  }
464 
465  // Set it, if it exists
466  if(fPHOSMatrix[mod])
467  {
468  if(fDebug > 1 )
469  fPHOSMatrix[mod]->Print();
470 
471  fPHOSGeo->SetMisalMatrix(fPHOSMatrix[mod],mod) ;
472  }
473  }// SM loop
474 
475  fPHOSGeoMatrixSet = kTRUE;//At least one, so good
476 
477  }//Load matrices
478  else if (!gGeoManager)
479  {
480  AliDebug(1,"Load PHOS misalignment matrices.");
481  if(!strcmp(inputEvent->GetName(),"AliESDEvent"))
482  {
483  for(Int_t mod = 0; mod < 5; mod++)
484  {
485  if( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod))
486  {
487  //printf("PHOS: mod %d, matrix %p\n",mod, ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod));
488  fPHOSGeo->SetMisalMatrix( ((AliESDEvent*)inputEvent)->GetPHOSMatrix(mod),mod) ;
489  }
490  } // loop over modules
491 
492  fPHOSGeoMatrixSet = kTRUE; //At least one so good
493  } // ESD as input
494  else
495  {
496  AliDebug(1,"Setting of EMCAL transformation matrixes for AODs not implemented yet. \n Import geometry.root file");
497  } // AOD as input
498  } // get matrix from data
499  else if(gGeoManager)
500  {
501  fPHOSGeoMatrixSet = kTRUE;
502  }
503  }//PHOS geo and geoManager was not set
504 }
505 
506 //______________________________________________________________________________________
509 //______________________________________________________________________________________
510 Bool_t AliCalorimeterUtils::AreNeighbours(Int_t calo, Int_t absId1, Int_t absId2 ) const
511 {
512  Bool_t areNeighbours = kFALSE ;
513 
514  Int_t iRCU1 = -1, irow1 = -1, icol1 = -1;
515  Int_t iRCU2 = -1, irow2 = -1, icol2 = -1;
516 
517  Int_t rowdiff = 0, coldiff = 0;
518 
519  Int_t nSupMod1 = GetModuleNumberCellIndexes(absId1, calo, icol1, irow1, iRCU1);
520  Int_t nSupMod2 = GetModuleNumberCellIndexes(absId2, calo, icol2, irow2, iRCU2);
521 
522  if(calo==kEMCAL && nSupMod1!=nSupMod2)
523  {
524  // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1
525  // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
526  if(nSupMod1%2) icol1+=AliEMCALGeoParams::fgkEMCALCols;
527  else icol2+=AliEMCALGeoParams::fgkEMCALCols;
528  }
529 
530  rowdiff = TMath::Abs( irow1 - irow2 ) ;
531  coldiff = TMath::Abs( icol1 - icol2 ) ;
532 
533  //if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
534  if ((coldiff + rowdiff == 1 ))
535  areNeighbours = kTRUE ;
536 
537  return areNeighbours;
538 }
539 
540 //_________________________________________________________________________________________
543 //_________________________________________________________________________________________
544 Bool_t AliCalorimeterUtils::IsClusterSharedByTwoSuperModules(const AliEMCALGeometry * geom,
545  AliVCluster* cluster)
546 {
547  Int_t iSupMod = -1;
548  Int_t iSM0 = -1;
549  Int_t iTower = -1;
550  Int_t iIphi = -1;
551  Int_t iIeta = -1;
552  Int_t iphi = -1;
553  Int_t ieta = -1;
554 
555  for(Int_t iDigit = 0; iDigit < cluster->GetNCells(); iDigit++)
556  {
557  // Get from the absid the supermodule, tower and eta/phi numbers
558  geom->GetCellIndex(cluster->GetCellAbsId(iDigit),iSupMod,iTower,iIphi,iIeta);
559  geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi,iIeta, iphi,ieta);
560 
561  // Check if there are cells of different SM
562  if (iDigit == 0 ) iSM0 = iSupMod;
563  else if(iSupMod != iSM0)
564  {
565  if(iSupMod > 11 && iSupMod < 18)
566  AliWarning(Form("Cluster shared in 2 DCal: SM%d, SM%d??",iSupMod,iSM0));
567 
568  return kTRUE;
569  }
570  }
571 
572  return kFALSE;
573 }
574 
575 //______________________________________________________________________________
578 //______________________________________________________________________________
579 Bool_t AliCalorimeterUtils::CheckCellFiducialRegion(AliVCluster* cluster,
580  AliVCaloCells* cells) const
581 {
582  //If the distance to the border is 0 or negative just exit accept all clusters
583 
584  if ( cells->GetType()==AliVCaloCells::kEMCALCell && fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder() <= 0 ) return kTRUE;
585 
586  if ( cells->GetType()==AliVCaloCells::kPHOSCell && fNCellsFromPHOSBorder <= 0 ) return kTRUE;
587 
588  Int_t absIdMax = -1;
589  Float_t ampMax = -1;
590 
591  for(Int_t i = 0; i < cluster->GetNCells() ; i++)
592  {
593  Int_t absId = cluster->GetCellAbsId(i) ;
594  Float_t amp = cells->GetCellAmplitude(absId);
595 
596  if(amp > ampMax)
597  {
598  ampMax = amp;
599  absIdMax = absId;
600  }
601  }
602 
603  AliDebug(1,Form("Cluster Max AbsId %d, Cell Energy %2.2f, Cluster Energy %2.2f",
604  absIdMax, ampMax, cluster->E()));
605 
606  if ( absIdMax == -1 ) return kFALSE;
607 
608  // Check if the cell is close to the borders:
609  Bool_t okrow = kFALSE;
610  Bool_t okcol = kFALSE;
611 
612  if ( cells->GetType() == AliVCaloCells::kEMCALCell )
613  {
614  // It should be the same as AliEMCALRecoUtils::CheckCellFiducialRegion()
615  // Why not calling it?
616 
617  Int_t iTower = -1, iIphi = -1, iIeta = -1, iphi = -1, ieta = -1, iSM = -1;
618 
619  fEMCALGeo->GetCellIndex(absIdMax,iSM,iTower,iIphi,iIeta);
620 
621  fEMCALGeo->GetCellPhiEtaIndexInSModule(iSM,iTower,iIphi, iIeta,iphi,ieta);
622 
623  if(iSM < 0 || iphi < 0 || ieta < 0 )
624  {
625  AliFatal(Form("Negative value for super module: %d, or cell ieta: %d, or cell iphi: %d, check EMCAL geometry name",iSM,ieta,iphi));
626  }
627 
628  // Check rows/phi
629  Int_t nborder = fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder();
630 
631  if ( iSM < 10 || (iSM > 11 && iSM < 18) ) // Full EMCal (SM0-9) and DCal 2/3 (SM12-17)
632  {
633  if(iphi >= nborder && iphi < 24-nborder) okrow = kTRUE;
634  }
635  else // 1/3 SMs (SM10-11, SM18-19)
636  {
637  if(iphi >= nborder && iphi < 8-nborder) okrow = kTRUE;
638  }
639 
640  // Check columns/eta
641 
642  // Remove all borders if IsEMCALNoBorderAtEta0 or DCal SMs(12-17)
643  if(!fEMCALRecoUtils->IsEMCALNoBorderAtEta0() || (iSM > 11 && iSM < 18))
644  {
645  if(ieta > nborder && ieta < 48-nborder) okcol =kTRUE;
646  }
647  else // Do not remove borders close at eta = 0 for Full EMCal SMs and 1/3 EMCal
648  {
649  if ( iSM%2 == 0 )
650  {
651  if(ieta >= nborder) okcol = kTRUE;
652  }
653  else
654  {
655  if(ieta < 48-nborder) okcol = kTRUE;
656  }
657  } // eta 0 not checked
658 
659  AliDebug(1,Form("EMCAL Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ? ok row %d, ok column %d",
660  nborder, ieta, iphi, iSM,okrow,okcol));
661 
662  }//EMCAL
663  else if ( cells->GetType() == AliVCaloCells::kPHOSCell )
664  {
665  Int_t relId[4];
666  Int_t irow = -1, icol = -1;
667  fPHOSGeo->AbsToRelNumbering(absIdMax,relId);
668 
669  if (relId[1] != 0 ) return kFALSE; // skip CPV only PHOS
670 
671  irow = relId[2];
672  icol = relId[3];
673  //imod = relId[0]-1;
674 
675  if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
676  if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
677 
678  AliDebug(1,Form("PHOS Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ? ok row %d, ok column %d",
679  fNCellsFromPHOSBorder, icol, irow, relId[0]-1,okrow,okcol));
680  }//PHOS
681 
682  if (okcol && okrow) return kTRUE;
683  else return kFALSE;
684 }
685 
686 //________________________________________________________________________________________________________
689 //________________________________________________________________________________________________________
690 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(Int_t calorimeter, UShort_t* cellList, Int_t nCells)
691 {
692  if (!fRemoveBadChannels) return kFALSE;
693 
694  //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
695  if(calorimeter == kEMCAL && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
696  if(calorimeter == kPHOS && !fPHOSBadChannelMap) return kFALSE;
697 
698  Int_t icol = -1;
699  Int_t irow = -1;
700  Int_t imod = -1;
701  for(Int_t iCell = 0; iCell<nCells; iCell++)
702  {
703  // Get the column and row
704  if ( calorimeter == kEMCAL )
705  {
706  return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
707  }
708  else if ( calorimeter == kPHOS )
709  {
710  Int_t relId[4];
711  fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
712 
713  if (relId[1] != 0 ) return kTRUE; // skip CPV only PHOS
714 
715  irow = relId[2];
716  icol = relId[3];
717  imod = relId[0]-1;
718 
719  if ( fPHOSBadChannelMap->GetEntries() <= imod ) continue;
720 
721  //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
722  if ( GetPHOSChannelStatus(imod, irow, icol) ) return kTRUE;
723  }
724  else return kFALSE;
725  } // cell cluster loop
726 
727  return kFALSE;
728 }
729 
730 //_______________________________________________________________
732 //_______________________________________________________________
734 {
735  clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
736 }
737 
738 //_______________________________________________________________
746 //______________________________________________________________________________
747 void AliCalorimeterUtils::GetEMCALSubregion(AliVCluster * clus, AliVCaloCells * cells,
748  Int_t & regEta, Int_t & regPhi) const
749 {
750  regEta = regPhi = -1 ;
751 
752  if(!clus->IsEMCAL()) return ;
753 
754  Int_t icol = -1, irow = -1, iRCU = -1;
755  Float_t clusterFraction = 0;
756 
757  Int_t absId = GetMaxEnergyCell(cells,clus,clusterFraction);
758 
759  Int_t sm = GetModuleNumberCellIndexes(absId,kEMCAL,icol,irow,iRCU);
760 
761  // Shift by 48 to for impair SM
762  if( sm%2 == 1) icol+=AliEMCALGeoParams::fgkEMCALCols;
763 
764  // Avoid borders
765  if(icol < 2 || icol > 93 || irow < 2 || irow > 21) return;
766 
767  //
768  // Eta regions
769  //
770 
771  // Region 0: center of SM ~0.18<|eta|<0.55
772  if ( icol > 9 && icol < 34 ) regEta = 0;
773  else if ( icol > 62 && icol < 87 ) regEta = 0;
774 
775  // Region 3: frames ~0.1<|eta|<~0.22 ~0.51<|eta|<0.62
776 
777  else if ( icol <= 9 && icol >= 5 ) regEta = 3;
778  else if ( icol <= 38 && icol >= 34 ) regEta = 3;
779  else if ( icol <= 62 && icol >= 58 ) regEta = 3;
780  else if ( icol <= 91 && icol >= 87 ) regEta = 3;
781 
782  // Region 1: |eta| < ~0.15
783 
784  else if ( icol < 58 && icol > 38 ) regEta = 1 ;
785 
786  // Region 2: |eta| > ~0.6
787 
788  else regEta = 2 ;
789 
790  //
791  // Phi regions
792  //
793 
794  if ( irow >= 2 && irow <= 5 ) regPhi = 0; // External
795  else if ( irow >= 18 && irow <= 21 ) regPhi = 0; // External
796  else if ( irow >= 6 && irow <= 9 ) regPhi = 1; // Mid
797  else if ( irow >= 14 && irow <= 17 ) regPhi = 1; // Mid
798  else regPhi = 2; //10-13 Central
799 
800 }
801 
802 //________________________________________________________________________________________
804 //________________________________________________________________________________________
805 Int_t AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells * cells,
806  AliVCluster * clu,
807  Float_t & clusterFraction) const
808 {
809  if( !clu || !cells )
810  {
811  AliInfo("Cluster or cells pointer is null!");
812  return -1;
813  }
814 
815  Double_t eMax =-1.;
816  Double_t eTot = 0.;
817  Double_t eCell =-1.;
818  Float_t fraction = 1.;
819  Float_t recalFactor = 1.;
820  Int_t cellAbsId =-1 , absId =-1 ;
821  Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
822 
823  Int_t calo = kEMCAL;
824  if(clu->IsPHOS()) calo = kPHOS ;
825 
826  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
827  {
828  cellAbsId = clu->GetCellAbsId(iDig);
829 
830  fraction = clu->GetCellAmplitudeFraction(iDig);
831  if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
832 
833  iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
834 
835  if(IsRecalibrationOn())
836  {
837  if(calo == kEMCAL) recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
838  else recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
839  }
840 
841  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
842 
843  if(eCell > eMax)
844  {
845  eMax = eCell;
846  absId = cellAbsId;
847  }
848 
849  eTot+=eCell;
850 
851  }// cell loop
852 
853  if(eTot > 0.1)
854  clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
855  else
856  clusterFraction =1.;
857 
858  return absId;
859 }
860 
861 //___________________________________________________________________________________
865 //___________________________________________________________________________________
866 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(AliVCluster* cluster,
867  AliVEvent* event, Int_t index) const
868 {
869  AliVTrack *track = 0x0;
870 
871  //
872  // EMCAL case only when matching is recalculated
873  //
874  if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
875  {
876  Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
877  //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
878 
879  if(trackIndex < 0 )
880  AliInfo(Form("Wrong track index %d, from recalculation", trackIndex));
881  else
882  track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
883 
884  return track ;
885  }
886 
887  //
888  // Normal case, get info from ESD or AOD
889  //
890 
891  // No tracks matched
892  if( cluster->GetNTracksMatched() < 1 ) return 0x0;
893 
894  // At least one match
895  Int_t iTrack = 0; // only one match for AODs with index 0.
896 
897  // ESDs
898  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
899  {
900  if( index >= 0 ) iTrack = index;
901  else iTrack = ((AliESDCaloCluster*)cluster)->GetTracksMatched()->At(0); //cluster->GetTrackMatchedIndex();
902 
903  track = dynamic_cast<AliVTrack*> ( event->GetTrack(iTrack) );
904  }
905  else // AODs
906  {
907  if( index > 0 ) iTrack = index;
908 
909  track = dynamic_cast<AliVTrack*>( cluster->GetTrackMatched(iTrack) );
910  }
911 
912  return track ;
913 }
914 
919 Float_t AliCalorimeterUtils::GetMCECellClusFracCorrection(Float_t eCell, Float_t eCluster) const
920 {
921  if( eCluster <= 0 || eCluster < eCell )
922  {
923  AliWarning(Form("Bad values eCell=%f, eCluster %f",eCell,eCluster));
924  return 1;
925  }
926 
927  Float_t frac = eCell / eCluster;
928 
929  Float_t correction = fMCECellClusFracCorrParam[0] +
930  TMath::Exp( frac*fMCECellClusFracCorrParam[2]+fMCECellClusFracCorrParam[1] ) +
931  fMCECellClusFracCorrParam[3]/TMath::Sqrt(frac);
932 
933 // printf("AliCalorimeterUtils::GetMCECellClusFracCorrection(eCell=%f, eCluster %f, frac %f) = %f\n",eCell, eCluster, frac, correction);
934 // printf("\t %2.2f + TMath::Exp( %2.3f*%2.2f + %2.2f ) + %2.2f/TMath::Sqrt(%2.3f)) = %f\n",
935 // fMCECellClusFracCorrParam[0],frac,fMCECellClusFracCorrParam[2],fMCECellClusFracCorrParam[1],fMCECellClusFracCorrParam[3], frac, correction);
936 
937  return correction;
938 }
939 
940 //_____________________________________________________________________________________________________
942 //_____________________________________________________________________________________________________
943 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
944 {
945  Int_t absId = -1;
946 
947  if(particle->GetDetectorTag()==kEMCAL)
948  {
949  fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
950 
951  AliDebug(2,Form("EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d",
952  particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId)));
953 
954  return fEMCALGeo->GetSuperModuleNumber(absId) ;
955  } // EMCAL
956  else if ( particle->GetDetectorTag() == kPHOS )
957  {
958  // In case we use the MC reader, the input are TParticles,
959  // in this case use the corresponing method in PHOS Geometry to get the particle.
960  if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
961  {
962  Int_t mod =-1;
963  Double_t z = 0., x=0.;
964  TParticle* primary = 0x0;
965  AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
966 
967  if(stack)
968  {
969  primary = stack->Particle(particle->GetCaloLabel(0));
970  }
971  else
972  {
973  AliFatal("Stack not available, stop!");
974  }
975 
976  if(primary)
977  {
978  fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
979  }
980  else
981  {
982  AliFatal("Primary not available, stop!");
983  }
984  return mod;
985  }
986  // Input are ESDs or AODs, get the PHOS module number like this.
987  else
988  {
989  //FIXME
990  //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
991  //return GetModuleNumber(cluster);
992  //MEFIX
993  return -1;
994  }
995  } // PHOS
996 
997  return -1;
998 }
999 
1000 //_____________________________________________________________________
1002 //_____________________________________________________________________
1003 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
1004 {
1005  if(!cluster)
1006  {
1007  AliDebug(1,"AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
1008 
1009  return -1;
1010  }
1011 
1012  if ( cluster->GetNCells() <= 0 ) return -1;
1013 
1014  Int_t absId = cluster->GetCellAbsId(0);
1015 
1016  if ( absId < 0 ) return -1;
1017 
1018  if( cluster->IsEMCAL() )
1019  {
1020  AliDebug(2,Form("EMCAL absid %d, SuperModule %d",absId, fEMCALGeo->GetSuperModuleNumber(absId)));
1021 
1022  return fEMCALGeo->GetSuperModuleNumber(absId) ;
1023  } // EMCAL
1024  else if ( cluster->IsPHOS() )
1025  {
1026  Int_t relId[4];
1027  fPHOSGeo->AbsToRelNumbering(absId,relId);
1028 
1029  if (relId[1] != 0 ) return -1; // skip CPV only PHOS
1030 
1031  AliDebug(2,Form("PHOS absid %d Module %d",absId, relId[0]-1));
1032 
1033  return relId[0]-1;
1034  } // PHOS
1035 
1036  return -1;
1037 }
1038 
1039 //___________________________________________________________________________________________________
1041 //___________________________________________________________________________________________________
1043  Int_t & icol, Int_t & irow, Int_t & iRCU) const
1044 {
1045  Int_t imod = -1;
1046 
1047  if ( absId < 0 ) return -1 ;
1048 
1049  if ( calo == kEMCAL )
1050  {
1051  Int_t iTower = -1, iIphi = -1, iIeta = -1;
1052  fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
1053  fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
1054 
1055  if(imod < 0 || irow < 0 || icol < 0 )
1056  {
1057  AliFatal(Form("Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name",imod,icol,irow));
1058  }
1059 
1060  // In case of DCal C side, shift columns to match offline/online numbering
1061  // when calculating the DDL for Run2
1062  // See AliEMCALRawUtils::Digits2Raw and Raw2Digits.
1063  Int_t ico2 = icol ;
1064  if ( imod == 13 || imod == 15 || imod == 17 ) ico2 += 16;
1065 
1066  // RCU / DDL
1067  if(imod < 10 || (imod > 11 && imod < 18)) // (EMCAL Full || DCAL 2/3)
1068  {
1069  // RCU0 / DDL0
1070  if ( 0 <= irow && irow < 8 ) iRCU = 0; // first cable row
1071  else if ( 8 <= irow && irow < 16 &&
1072  0 <= ico2 && ico2 < 24 ) iRCU = 0; // first half;
1073  //second cable row
1074 
1075  // RCU1 / DDL1
1076  else if ( 8 <= irow && irow < 16 &&
1077  24 <= ico2 && ico2 < 48 ) iRCU = 1; // second half;
1078  //second cable row
1079  else if ( 16 <= irow && irow < 24 ) iRCU = 1; // third cable row
1080 
1081  if ( imod%2 == 1 ) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1082  }
1083  else
1084  {
1085  // 1/3 SM have one single SRU, just assign RCU/DDL 0
1086  iRCU = 0 ;
1087  }
1088 
1089  if ( iRCU < 0 )
1090  AliFatal(Form("Wrong EMCAL RCU number = %d", iRCU));
1091 
1092  return imod ;
1093  } // EMCAL
1094  else // PHOS
1095  {
1096  Int_t relId[4];
1097  fPHOSGeo->AbsToRelNumbering(absId,relId);
1098 
1099  if (relId[1] != 0 ) return -1; // skip CPV only PHOS
1100 
1101  irow = relId[2];
1102  icol = relId[3];
1103  imod = relId[0]-1;
1104  iRCU= (Int_t)(relId[2]-1)/16 ;
1105 
1106  //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
1107 
1108  if ( iRCU >= 4 )
1109  AliFatal(Form("Wrong PHOS RCU number = %d", iRCU));
1110 
1111  return imod;
1112  } // PHOS
1113 
1114  return -1;
1115 }
1116 
1117 //___________________________________________________________________________________________
1119 //___________________________________________________________________________________________
1120 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1121 {
1122  const Int_t nc = cluster->GetNCells();
1123 
1124  Int_t absIdList[nc];
1125  Float_t maxEList[nc];
1126 
1127  Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1128 
1129  return nMax;
1130 }
1131 
1132 //___________________________________________________________________________________________
1134 //___________________________________________________________________________________________
1135 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1136  Int_t *absIdList, Float_t *maxEList)
1137 {
1138  Int_t iDigitN = 0 ;
1139  Int_t iDigit = 0 ;
1140  Int_t absId1 = -1 ;
1141  Int_t absId2 = -1 ;
1142  const Int_t nCells = cluster->GetNCells();
1143 
1144  Float_t eCluster = RecalibrateClusterEnergy(cluster, cells);// recalculate cluster energy, avoid non lin correction.
1145 
1146  Float_t simuTotWeight = 0;
1148  {
1149  simuTotWeight = RecalibrateClusterEnergyWeightCell(cluster, cells,eCluster);// same but apply a weight
1150  simuTotWeight/= eCluster;
1151  }
1152 
1153  Int_t calorimeter = kEMCAL;
1154  if(!cluster->IsEMCAL()) calorimeter = kPHOS;
1155 
1156  //printf("cluster : ncells %d \n",nCells);
1157 
1158  Float_t emax = 0;
1159  Int_t idmax =-1;
1160  for(iDigit = 0; iDigit < nCells ; iDigit++)
1161  {
1162  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
1163  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1164  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1165 
1167  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1168 
1169  if( en > emax )
1170  {
1171  emax = en ;
1172  idmax = absIdList[iDigit] ;
1173  }
1174  //Int_t icol = -1, irow = -1, iRCU = -1;
1175  //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1176  //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
1177  }
1178 
1179  for(iDigit = 0 ; iDigit < nCells; iDigit++)
1180  {
1181  if( absIdList[iDigit] >= 0 )
1182  {
1183  absId1 = cluster->GetCellsAbsId()[iDigit];
1184 
1185  Float_t en1 = cells->GetCellAmplitude(absId1);
1186  RecalibrateCellAmplitude(en1,calorimeter,absId1);
1187 
1189  en1*=GetMCECellClusFracCorrection(en1,eCluster)/simuTotWeight;
1190 
1191  //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
1192 
1193  for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1194  {
1195  absId2 = cluster->GetCellsAbsId()[iDigitN] ;
1196 
1197  if(absId2==-1 || absId2==absId1) continue;
1198 
1199  //printf("\t %d : absIDj %d\n",iDigitN, absId2);
1200 
1201  Float_t en2 = cells->GetCellAmplitude(absId2);
1202  RecalibrateCellAmplitude(en2,calorimeter,absId2);
1203 
1205  en2*=GetMCECellClusFracCorrection(en2,eCluster)/simuTotWeight;
1206 
1207  //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
1208 
1209  if ( AreNeighbours(calorimeter, absId1, absId2) )
1210  {
1211  // printf("\t \t Neighbours \n");
1212  if ( en1 > en2 )
1213  {
1214  absIdList[iDigitN] = -1 ;
1215  //printf("\t \t indexN %d not local max\n",iDigitN);
1216  // but may be digit too is not local max ?
1217  if(en1 < en2 + fLocMaxCutEDiff) {
1218  //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
1219  absIdList[iDigit] = -1 ;
1220  }
1221  }
1222  else
1223  {
1224  absIdList[iDigit] = -1 ;
1225  //printf("\t \t index %d not local max\n",iDigitN);
1226  // but may be digitN too is not local max ?
1227  if(en1 > en2 - fLocMaxCutEDiff)
1228  {
1229  absIdList[iDigitN] = -1 ;
1230  //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
1231  }
1232  }
1233  } // if Are neighbours
1234  //else printf("\t \t NOT Neighbours \n");
1235  } // while digitN
1236  } // slot not empty
1237  } // while digit
1238 
1239  iDigitN = 0 ;
1240  for(iDigit = 0; iDigit < nCells; iDigit++)
1241  {
1242  if( absIdList[iDigit] >= 0 )
1243  {
1244  absIdList[iDigitN] = absIdList[iDigit] ;
1245 
1246  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1247  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1248 
1250  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1251 
1252  if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1253 
1254  maxEList[iDigitN] = en ;
1255 
1256  //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1257  iDigitN++ ;
1258  }
1259  }
1260 
1261  if ( iDigitN == 0 )
1262  {
1263  AliDebug(1,Form("No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f",
1264  idmax,emax,cluster->E()));
1265  iDigitN = 1 ;
1266  maxEList[0] = emax ;
1267  absIdList[0] = idmax ;
1268  }
1269 
1270 
1271  AliDebug(1,Form("In cluster E %2.2f (wth non lin. %2.2f), M02 %2.2f, M20 %2.2f, N maxima %d",
1272  cluster->E(),eCluster, cluster->GetM02(),cluster->GetM20(), iDigitN));
1273 
1274 // if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1275 // {
1276 // printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1277 // }
1278 
1279  return iDigitN ;
1280 }
1281 
1282 //____________________________________
1284 //____________________________________
1286 {
1287  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1288  {
1289  AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null");
1290  return TString("");
1291  }
1292 
1293  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1294  {
1295  AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null");
1296  return TString("");
1297  }
1298 
1299  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1300  if (pass.Contains("ass1")) return TString("pass1");
1301  else if (pass.Contains("ass2")) return TString("pass2");
1302  else if (pass.Contains("ass3")) return TString("pass3");
1303  else if (pass.Contains("ass4")) return TString("pass4");
1304  else if (pass.Contains("ass5")) return TString("pass5");
1305  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1306  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1307  {
1308  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1309  return TString("pass1");
1310  }
1311 
1312  // No condition fullfilled, give a default value
1313  AliInfo("Pass number string not found");
1314  return TString("");
1315 }
1316 
1317 //________________________________________
1319 //________________________________________
1321 {
1322  fEMCALGeoName = "";
1323  fPHOSGeoName = "";
1324 
1325  fEMCALGeoMatrixSet = kFALSE;
1326  fPHOSGeoMatrixSet = kFALSE;
1327 
1328  fRemoveBadChannels = kFALSE;
1329 
1331 
1332  fLocMaxCutE = 0.1 ;
1333  fLocMaxCutEDiff = 0.0 ;
1334 
1335  // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1336  // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1337  // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1338  // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1339  // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1340  // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
1341 
1342  fOADBSet = kFALSE;
1343  fOADBForEMCAL = kTRUE ;
1344  fOADBForPHOS = kFALSE;
1345 
1346  fOADBFilePathEMCAL = "$ALICE_PHYSICS/OADB/EMCAL" ;
1347  fOADBFilePathPHOS = "$ALICE_PHYSICS/OADB/PHOS" ;
1348 
1349  fImportGeometryFromFile = kTRUE;
1351 
1352  fNSuperModulesUsed = 22;
1353 
1354  fMCECellClusFracCorrParam[0] = 0.78;
1355  fMCECellClusFracCorrParam[1] =-1.8;
1356  fMCECellClusFracCorrParam[2] =-6.3;
1357  fMCECellClusFracCorrParam[3] = 0.014;
1358 }
1359 
1360 //_____________________________________________________
1362 //_____________________________________________________
1364 {
1365  AliDebug(1,"Init bad channel map");
1366 
1367  // In order to avoid rewriting the same histograms
1368  Bool_t oldStatus = TH1::AddDirectoryStatus();
1369  TH1::AddDirectory(kFALSE);
1370 
1371  fPHOSBadChannelMap = new TObjArray(5);
1372  for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),
1373  Form("PHOS_BadMap_mod%d",i),
1374  64, 0, 64, 56, 0, 56));
1375 
1376  fPHOSBadChannelMap->SetOwner(kTRUE);
1377  fPHOSBadChannelMap->Compress();
1378 
1379  //In order to avoid rewriting the same histograms
1380  TH1::AddDirectory(oldStatus);
1381 }
1382 
1383 //______________________________________________________
1385 //______________________________________________________
1387 {
1388  AliDebug(1,"Init recalibration map");
1389 
1390  // In order to avoid rewriting the same histograms
1391  Bool_t oldStatus = TH1::AddDirectoryStatus();
1392  TH1::AddDirectory(kFALSE);
1393 
1394  fPHOSRecalibrationFactors = new TObjArray(5);
1395  for (int i = 0; i < 5; i++)
1396  {
1397  fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),
1398  Form("PHOSRecalFactors_Mod%d",i),
1399  64, 0, 64, 56, 0, 56));
1400  }
1401 
1402  // Init the histograms with 1
1403  for (Int_t m = 0; m < 5; m++)
1404  {
1405  for (Int_t i = 0; i < 56; i++)
1406  {
1407  for (Int_t j = 0; j < 64; j++)
1408  {
1410  }
1411  }
1412  }
1413 
1414  fPHOSRecalibrationFactors->SetOwner(kTRUE);
1415  fPHOSRecalibrationFactors->Compress();
1416 
1417  // In order to avoid rewriting the same histograms
1418  TH1::AddDirectory(oldStatus);
1419 }
1420 
1425 {
1426  if (fEMCALGeo) return;
1427 
1428  AliDebug(1,Form(" for run=%d",fRunNumber));
1429 
1430  if(fEMCALGeoName=="")
1431  {
1432  fEMCALGeo = AliEMCALGeometry::GetInstanceFromRunNumber(fRunNumber);
1433  AliInfo(Form("Get EMCAL geometry name to <%s> for run %d",fEMCALGeo->GetName(),fRunNumber));
1434  }
1435  else
1436  {
1437  fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
1438  AliInfo(Form("Set EMCAL geometry name to <%s>",fEMCALGeoName.Data()));
1439  }
1440 
1441  // Init geometry, I do not like much to do it like this ...
1442  if(fImportGeometryFromFile && !gGeoManager)
1443  {
1444  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1445  {
1446  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1447  if (fRunNumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1448  else if(fRunNumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1449  else if(fRunNumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1450  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >= 2015
1451  }
1452 
1453  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1454 
1455  TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1456  }
1457  else if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1458 }
1459 
1464 {
1465  if (fPHOSGeo) return;
1466 
1467  AliDebug(1,Form(" for run=%d",fRunNumber));
1468 
1469  if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1470 
1471  fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
1472 
1473  //if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1474 }
1475 
1476 //______________________________________________________________________________________________
1477 // Check that a MC ESD is in the calorimeter acceptance
1478 //______________________________________________________________________________________________
1479 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(Int_t calo, TParticle* particle)
1480 {
1481  if(!particle || (calo!=kEMCAL && calo!=kPHOS)) return kFALSE ;
1482 
1483  if( (!IsPHOSGeoMatrixSet () && calo == kPHOS ) ||
1484  (!IsEMCALGeoMatrixSet() && calo == kEMCAL) )
1485  {
1486  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1487  return kFALSE ;
1488  }
1489 
1490  if(calo == kPHOS )
1491  {
1492  Int_t mod = 0 ;
1493  Double_t x = 0, z = 0 ;
1494  return GetPHOSGeometry()->ImpactOnEmc( particle, mod, z, x);
1495  }
1496  else if(calo == kEMCAL)
1497  {
1498  Int_t absID = 0 ;
1499  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(),absID);
1500  if(ok)
1501  {
1502  Int_t icol = -1, irow = -1, iRCU = -1;
1503  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1504  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1505  if(status > 0) ok = kFALSE;
1506  }
1507 
1508  return ok ;
1509  }
1510 
1511  return kFALSE ;
1512 }
1513 
1514 //______________________________________________________________________________________________________
1516 //______________________________________________________________________________________________________
1517 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(Int_t calo, AliAODMCParticle* particle)
1518 {
1519  if(!particle || (calo!=kEMCAL && calo!=kPHOS)) return kFALSE ;
1520 
1521  if( (!IsPHOSGeoMatrixSet () && calo == kPHOS ) ||
1522  (!IsEMCALGeoMatrixSet() && calo == kEMCAL) )
1523  {
1524  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1525  return kFALSE ;
1526  }
1527 
1528  Float_t phi = particle->Phi();
1529  if(phi < 0) phi+=TMath::TwoPi();
1530 
1531  if(calo == kPHOS )
1532  {
1533  Int_t mod = 0 ;
1534  Double_t x = 0, z = 0 ;
1535  Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1536  return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), phi, mod, z, x) ;
1537  }
1538  else if(calo == kEMCAL)
1539  {
1540  Int_t absID = 0 ;
1541  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),phi,absID);
1542  if(ok)
1543  {
1544  Int_t icol = -1, irow = -1, iRCU = -1;
1545  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1546  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1547  if(status > 0) ok = kFALSE;
1548  }
1549 
1550  return ok ;
1551  }
1552 
1553  return kFALSE ;
1554 }
1555 
1556 //_____________________________________________________________________________________________________
1557 // Check that a TLorentzVector is in the calorimeter acceptance, give the cell number where it hit.
1558 //_____________________________________________________________________________________________________
1559 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(Int_t calo, Float_t eta, Float_t theta,
1560  Float_t phiOrg, Int_t & absID)
1561 {
1562  if(calo!=kEMCAL && calo!=kPHOS) return kFALSE ;
1563 
1564  if( (!IsPHOSGeoMatrixSet () && calo == kPHOS ) ||
1565  (!IsEMCALGeoMatrixSet() && calo == kEMCAL) )
1566  {
1567  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1568  return kFALSE ;
1569  }
1570 
1571  Float_t phi = phiOrg;
1572  if(phi < 0) phi+=TMath::TwoPi();
1573 
1574  if(calo == kPHOS )
1575  {
1576  Int_t mod = 0 ;
1577  Double_t x = 0, z = 0 ;
1578  Double_t vtx[]={0,0,0} ;
1579  return GetPHOSGeometry()->ImpactOnEmc(vtx, theta, phi, mod, z, x) ;
1580  }
1581  else if(calo == kEMCAL)
1582  {
1583  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta,phi,absID);
1584  if(ok)
1585  {
1586  Int_t icol = -1, irow = -1, iRCU = -1;
1587  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1588  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1589  if(status > 0) ok = kFALSE;
1590  }
1591 
1592  return ok ;
1593  }
1594 
1595  return kFALSE ;
1596 }
1597 
1598 //_______________________________________________________________________
1601 //_______________________________________________________________________
1602 Bool_t AliCalorimeterUtils::MaskFrameCluster(Int_t iSM, Int_t ieta) const
1603 {
1604  Int_t icol = ieta;
1605  if ( iSM%2 ) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1606 
1608  {
1609  for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1610  {
1611  if(icol==fMaskCellColumns[imask]) return kTRUE;
1612  }
1613  }
1614 
1615  return kFALSE;
1616 }
1617 
1618 //_________________________________________________________
1620 //_________________________________________________________
1621 void AliCalorimeterUtils::Print(const Option_t * opt) const
1622 {
1623  if(! opt)
1624  return;
1625 
1626  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1627  printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1628  printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
1629  fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
1630  if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
1631  printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
1632  printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1633  printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
1634  printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1635 
1636  printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1637  printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1638 
1639  printf(" \n") ;
1640 }
1641 
1642 //_____________________________________________________________________________________________
1644 //_____________________________________________________________________________________________
1645 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp, Int_t calo, Int_t id) const
1646 {
1647  Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1648  Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1649 
1650  if (IsRecalibrationOn())
1651  {
1652  if(calo == kPHOS)
1653  {
1654  amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1655  }
1656  else
1657  {
1658  amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1659  }
1660  }
1661 }
1662 
1663 //____________________________________________________________________________________________________
1665 //____________________________________________________________________________________________________
1666 void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, Int_t calo, Int_t id, Int_t bc) const
1667 {
1668  if ( calo == kEMCAL && GetEMCALRecoUtils()->IsTimeRecalibrationOn() )
1669  {
1670  GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1671  }
1672 }
1673 
1674 //__________________________________________________________________________
1676 //__________________________________________________________________________
1677 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1678  AliVCaloCells * cells)
1679 {
1680  // Initialize some used variables
1681  Float_t frac = 0., energy = 0.;
1682 
1683  if(cells)
1684  {
1685  //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1686 
1687  UShort_t * index = cluster->GetCellsAbsId() ;
1688  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1689 
1690  Int_t ncells = cluster->GetNCells();
1691 
1692  Int_t calo = kEMCAL;
1693  if(cluster->IsPHOS()) calo = kPHOS ;
1694 
1695  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1696  for(Int_t icell = 0; icell < ncells; icell++)
1697  {
1698  Int_t absId = index[icell];
1699 
1700  frac = fraction[icell];
1701  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1702 
1703  Float_t amp = cells->GetCellAmplitude(absId);
1704  RecalibrateCellAmplitude(amp,calo, absId);
1705 
1706  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy: before cal %f; after cal %f",
1707  calo,frac,cells->GetCellAmplitude(absId),amp));
1708 
1709  energy += amp*frac;
1710  }
1711 
1712  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
1713 
1714  } // cells available
1715  else
1716  {
1717  AliFatal("Cells pointer does not exist!");
1718  }
1719 
1720  return energy;
1721 }
1722 
1723 //_______________________________________________________________________________________________________
1726 //_______________________________________________________________________________________________________
1728  AliVCaloCells * cells, Float_t energyOrg)
1729 {
1730  //Initialize some used variables
1731  Float_t frac = 0., energy = 0.;
1732 
1733  if(cells)
1734  {
1735  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1736 
1737  UShort_t * index = cluster->GetCellsAbsId() ;
1738  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1739 
1740  Int_t ncells = cluster->GetNCells();
1741 
1742  Int_t calo = kEMCAL;
1743  if(cluster->IsPHOS()) calo = kPHOS ;
1744 
1745  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1746  for(Int_t icell = 0; icell < ncells; icell++)
1747  {
1748  Int_t absId = index[icell];
1749 
1750  frac = fraction[icell];
1751  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1752 
1753  Float_t amp = cells->GetCellAmplitude(absId);
1754  RecalibrateCellAmplitude(amp,calo, absId);
1755 
1756  amp*=GetMCECellClusFracCorrection(amp,energyOrg);
1757 
1758  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy %f",
1759  calo,frac,cells->GetCellAmplitude(absId)));
1760 
1761  energy += amp*frac;
1762  }
1763 
1764  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
1765  } // cells available
1766  else
1767  {
1768  AliFatal("Cells pointer does not exist!");
1769  }
1770 
1771  return energy;
1772 }
1773 
1774 //__________________________________________________________________________________________
1777 //__________________________________________________________________________________________
1778 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1779 {
1780  fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
1781 }
1782 
1783 //________________________________________________________________________________
1789 //________________________________________________________________________________
1791  TObjArray* clusterArray)
1792 {
1793  if (!fRecalculateMatching) return ;
1794 
1795  fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1796 
1797  Float_t dZ = 2000;
1798  Float_t dR = 2000;
1799 
1800  Int_t nClusters = event->GetNumberOfCaloClusters();
1801  if(clusterArray) nClusters = clusterArray->GetEntriesFast();
1802 
1803  AliVCluster * clus = 0;
1804 
1805  for (Int_t iclus = 0; iclus < nClusters ; iclus++)
1806  {
1807  if ( clusterArray ) clus = (AliVCluster*) clusterArray->At(iclus) ;
1808  else clus = event->GetCaloCluster(iclus) ;
1809 
1810  if (!clus->IsEMCAL()) continue ;
1811 
1812  //
1813  // Put track residuals in cluster
1814  //
1815  fEMCALRecoUtils->GetMatchedResiduals(clus->GetID(),dZ,dR);
1816 
1817  if ( TMath::Abs(clus->GetTrackDx()) < 500 )
1818  AliDebug(2,Form("Residuals (Old, New): z (%2.4f,%2.4f), x (%2.4f,%2.4f)\n",
1819  clus->GetTrackDz(),dZ,clus->GetTrackDx(),dR));
1820 
1821  clus->SetTrackDistance(dR,dZ);
1822 
1823  //
1824  // Remove old matches in cluster
1825  //
1826  if(clus->GetNTracksMatched() > 0)
1827  {
1828  if(!strcmp("AliESDCaloCluster",Form("%s",clus->ClassName())))
1829  {
1830  TArrayI arrayTrackMatched(0);
1831  ((AliESDCaloCluster*)clus)->AddTracksMatched(arrayTrackMatched);
1832  }
1833  else
1834  {
1835  for(Int_t iTrack = 0; iTrack < clus->GetNTracksMatched(); iTrack++)
1836  {
1837  ((AliAODCaloCluster*)clus)->RemoveTrackMatched((TObject*)((AliAODCaloCluster*)clus)->GetTrackMatched(iTrack));
1838  }
1839  }
1840  }
1841 
1842  //
1843  // Now put first track index in cluster.
1844  //
1845  Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(iclus);
1846  if ( trackIndex >= 0 )
1847  {
1848  if(!strcmp("AliESDCaloCluster",Form("%s",clus->ClassName())))
1849  {
1850  TArrayI arrayTrackMatched(1);
1851  arrayTrackMatched[0] = trackIndex;
1852  ((AliESDCaloCluster*)clus)->AddTracksMatched(arrayTrackMatched);
1853  }
1854  else
1855  {
1856  ((AliAODCaloCluster*)clus)->AddTrackMatched((TObject*)event->GetTrack(trackIndex));
1857  }
1858  }
1859 
1860  } // cluster loop
1861 }
1862 
1863 //___________________________________________________________________________
1876 //___________________________________________________________________________
1877 void AliCalorimeterUtils::SplitEnergy(Int_t absId1, Int_t absId2,
1878  AliVCluster* cluster,
1879  AliVCaloCells* cells,
1880  //Float_t & e1, Float_t & e2,
1881  AliAODCaloCluster* cluster1,
1882  AliAODCaloCluster* cluster2,
1883  Int_t nMax, Int_t eventNumber)
1884 {
1885  TH2F* hClusterMap = 0 ;
1886  TH2F* hClusterLocMax = 0 ;
1887  TH2F* hCluster1 = 0 ;
1888  TH2F* hCluster2 = 0 ;
1889 
1890  if(fPlotCluster)
1891  {
1892  hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1893  hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1894  hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1895  hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1896  hClusterMap ->SetXTitle("column");
1897  hClusterMap ->SetYTitle("row");
1898  hClusterLocMax ->SetXTitle("column");
1899  hClusterLocMax ->SetYTitle("row");
1900  hCluster1 ->SetXTitle("column");
1901  hCluster1 ->SetYTitle("row");
1902  hCluster2 ->SetXTitle("column");
1903  hCluster2 ->SetYTitle("row");
1904  }
1905 
1906  Int_t calorimeter = kEMCAL;
1907  if(cluster->IsPHOS())
1908  {
1909  calorimeter=kPHOS;
1910  AliWarning("Not supported for PHOS yet");
1911  return;
1912  }
1913 
1914  const Int_t ncells = cluster->GetNCells();
1915  Int_t absIdList[ncells];
1916 
1917  Float_t e1 = 0, e2 = 0 ;
1918  Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1919  Float_t eCluster = 0;
1920  Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1921  for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1922  {
1923  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1924 
1925  Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1926  RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1927  eCluster+=ec;
1928 
1929  if(fPlotCluster)
1930  {
1931  //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1932  sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1933  if(sm > -1 && sm < 12) // just to avoid compilation warning
1934  {
1935  if(icol > maxCol) maxCol = icol;
1936  if(icol < minCol) minCol = icol;
1937  if(irow > maxRow) maxRow = irow;
1938  if(irow < minRow) minRow = irow;
1939  hClusterMap->Fill(icol,irow,ec);
1940  }
1941  }
1942  }
1943 
1944  // Init counters and variables
1945  Int_t ncells1 = 1 ;
1946  UShort_t absIdList1[9] ;
1947  Double_t fracList1 [9] ;
1948  absIdList1[0] = absId1 ;
1949  fracList1 [0] = 1. ;
1950 
1951  Float_t ecell1 = cells->GetCellAmplitude(absId1);
1952  RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1953  e1 = ecell1;
1954 
1955  Int_t ncells2 = 1 ;
1956  UShort_t absIdList2[9] ;
1957  Double_t fracList2 [9] ;
1958  absIdList2[0] = absId2 ;
1959  fracList2 [0] = 1. ;
1960 
1961  Float_t ecell2 = cells->GetCellAmplitude(absId2);
1962  RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1963  e2 = ecell2;
1964 
1965  if(fPlotCluster)
1966  {
1967  Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1968  sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1969  hClusterLocMax->Fill(icol1,irow1,ecell1);
1970  sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1971  hClusterLocMax->Fill(icol2,irow2,ecell2);
1972  }
1973 
1974  // Very rough way to share the cluster energy
1975  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1976  Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1977  Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1978 
1979  for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1980  {
1981  Int_t absId = absIdList[iDigit];
1982 
1983  if ( absId==absId1 || absId==absId2 || absId < 0 ) continue;
1984 
1985  Float_t ecell = cells->GetCellAmplitude(absId);
1986  RecalibrateCellAmplitude(ecell, calorimeter, absId);
1987 
1988  if(AreNeighbours(calorimeter, absId1,absId ))
1989  {
1990  absIdList1[ncells1]= absId;
1991 
1992  if(AreNeighbours(calorimeter, absId2,absId ))
1993  {
1994  fracList1[ncells1] = shareFraction1;
1995  e1 += ecell*shareFraction1;
1996  }
1997  else
1998  {
1999  fracList1[ncells1] = 1.;
2000  e1 += ecell;
2001  }
2002 
2003  ncells1++;
2004 
2005  } // neigbour to cell1
2006 
2007  if(AreNeighbours(calorimeter, absId2,absId ))
2008  {
2009  absIdList2[ncells2]= absId;
2010 
2011  if(AreNeighbours(calorimeter, absId1,absId ))
2012  {
2013  fracList2[ncells2] = shareFraction2;
2014  e2 += ecell*shareFraction2;
2015  }
2016  else
2017  {
2018  fracList2[ncells2] = 1.;
2019  e2 += ecell;
2020  }
2021 
2022  ncells2++;
2023 
2024  } // neigbour to cell2
2025  }
2026 
2027  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",
2028  nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2));
2029 
2030  cluster1->SetE(e1);
2031  cluster2->SetE(e2);
2032 
2033  cluster1->SetNCells(ncells1);
2034  cluster2->SetNCells(ncells2);
2035 
2036  cluster1->SetCellsAbsId(absIdList1);
2037  cluster2->SetCellsAbsId(absIdList2);
2038 
2039  cluster1->SetCellsAmplitudeFraction(fracList1);
2040  cluster2->SetCellsAmplitudeFraction(fracList2);
2041 
2042  // Correct linearity
2043  CorrectClusterEnergy(cluster1) ;
2044  CorrectClusterEnergy(cluster2) ;
2045 
2046  if(calorimeter==kEMCAL)
2047  {
2048  GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
2049  GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
2050  }
2051 
2052  if(fPlotCluster)
2053  {
2054  //printf("Cells of cluster1: ");
2055  for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
2056  {
2057  //printf(" %d ",absIdList1[iDigit]);
2058 
2059  sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
2060 
2061  Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
2062  RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
2063 
2064  if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
2065  hCluster1->Fill(icol,irow,ecell*shareFraction1);
2066  else
2067  hCluster1->Fill(icol,irow,ecell);
2068  }
2069 
2070  //printf(" \n ");
2071  //printf("Cells of cluster2: ");
2072 
2073  for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
2074  {
2075  //printf(" %d ",absIdList2[iDigit]);
2076 
2077  sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
2078 
2079  Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
2080  RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
2081 
2082  if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) && absId2!=absIdList2[iDigit])
2083  hCluster2->Fill(icol,irow,ecell*shareFraction2);
2084  else
2085  hCluster2->Fill(icol,irow,ecell);
2086  }
2087  //printf(" \n ");
2088 
2089  gStyle->SetPadRightMargin(0.1);
2090  gStyle->SetPadLeftMargin(0.1);
2091  gStyle->SetOptStat(0);
2092  gStyle->SetOptFit(000000);
2093 
2094  if(maxCol-minCol > maxRow-minRow)
2095  {
2096  maxRow+= maxCol-minCol;
2097  }
2098  else
2099  {
2100  maxCol+= maxRow-minRow;
2101  }
2102 
2103  TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
2104  c->Divide(2,2);
2105  c->cd(1);
2106  gPad->SetGridy();
2107  gPad->SetGridx();
2108  gPad->SetLogz();
2109  hClusterMap ->SetAxisRange(minCol, maxCol,"X");
2110  hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
2111  hClusterMap ->Draw("colz TEXT");
2112  c->cd(2);
2113  gPad->SetGridy();
2114  gPad->SetGridx();
2115  gPad->SetLogz();
2116  hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
2117  hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
2118  hClusterLocMax ->Draw("colz TEXT");
2119  c->cd(3);
2120  gPad->SetGridy();
2121  gPad->SetGridx();
2122  gPad->SetLogz();
2123  hCluster1 ->SetAxisRange(minCol, maxCol,"X");
2124  hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
2125  hCluster1 ->Draw("colz TEXT");
2126  c->cd(4);
2127  gPad->SetGridy();
2128  gPad->SetGridx();
2129  gPad->SetLogz();
2130  hCluster2 ->SetAxisRange(minCol, maxCol,"X");
2131  hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
2132  hCluster2 ->Draw("colz TEXT");
2133 
2134  if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
2135  eventNumber,cluster->E(),nMax,ncells1,ncells2));
2136 
2137  delete c;
2138  delete hClusterMap;
2139  delete hClusterLocMax;
2140  delete hCluster1;
2141  delete hCluster2;
2142  }
2143 }
2144 
Bool_t fRecalculatePosition
Recalculate cluster position.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
AliEMCALGeometry * fEMCALGeo
! EMCAL geometry pointer.
void InitPHOSRecalibrationFactors()
Init PHOS recalibration factors.
Bool_t IsRecalculationOfClusterTrackMatchingOn() const
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)
AliPHOSGeoUtils * fPHOSGeo
! PHOS geometry pointer.
AliEMCALRecoUtils * GetEMCALRecoUtils() const
TString fileName
void SwitchOnDistToBadChannelRecalculation()
AliPHOSGeoUtils * GetPHOSGeometry() const
Bool_t fOADBForEMCAL
Get calibration from OADB for EMCAL.
Float_t GetPHOSChannelRecalibrationFactor(Int_t imod, Int_t iCol, Int_t iRow) const
void RecalculateClusterTrackMatching(AliVEvent *event, TObjArray *clusterArray=0x0)
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.
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:35
Int_t fDebug
Debugging level.
Bool_t IsRecalibrationOn() const
Bool_t IsMCParticleInCalorimeterAcceptance(Int_t calo, TParticle *particle)
virtual void Print(const Option_t *opt) const
Print some relevant parameters set for the analysis.
Int_t GetModuleNumber(AliAODPWG4Particle *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
Bool_t fRecalculateMatching
Recalculate cluster position.
Bool_t MaskFrameCluster(Int_t iSM, Int_t ieta) const
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.
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.
energy
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
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.
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.
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
Int_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow) const
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)
void SplitEnergy(Int_t absId1, Int_t absId2, AliVCluster *cluster, AliVCaloCells *cells, AliAODCaloCluster *cluster1, AliAODCaloCluster *cluster2, Int_t nMax, Int_t eventNumber=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...
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