AliPhysics  vAN-20151012 (2287573)
 All Classes Namespaces Files Functions Variables 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  irow = relId[2];
669  icol = relId[3];
670  //imod = relId[0]-1;
671 
672  if(irow >= fNCellsFromPHOSBorder && irow < 64-fNCellsFromPHOSBorder) okrow =kTRUE;
673  if(icol >= fNCellsFromPHOSBorder && icol < 56-fNCellsFromPHOSBorder) okcol =kTRUE;
674 
675  AliDebug(1,Form("PHOS Cluster in %d cells fiducial volume: ieta %d, iphi %d, SM %d ? ok row %d, ok column %d",
676  fNCellsFromPHOSBorder, icol, irow, relId[0]-1,okrow,okcol));
677  }//PHOS
678 
679  if (okcol && okrow) return kTRUE;
680  else return kFALSE;
681 }
682 
683 //________________________________________________________________________________________________________
686 //________________________________________________________________________________________________________
687 Bool_t AliCalorimeterUtils::ClusterContainsBadChannel(Int_t calorimeter, UShort_t* cellList, Int_t nCells)
688 {
689  if (!fRemoveBadChannels) return kFALSE;
690 
691  //printf("fEMCALBadChannelMap %p, fPHOSBadChannelMap %p \n",fEMCALBadChannelMap,fPHOSBadChannelMap);
692  if(calorimeter == kEMCAL && !fEMCALRecoUtils->GetEMCALChannelStatusMap(0)) return kFALSE;
693  if(calorimeter == kPHOS && !fPHOSBadChannelMap) return kFALSE;
694 
695  Int_t icol = -1;
696  Int_t irow = -1;
697  Int_t imod = -1;
698  for(Int_t iCell = 0; iCell<nCells; iCell++)
699  {
700  // Get the column and row
701  if ( calorimeter == kEMCAL )
702  {
703  return fEMCALRecoUtils->ClusterContainsBadChannel((AliEMCALGeometry*)fEMCALGeo,cellList,nCells);
704  }
705  else if ( calorimeter == kPHOS )
706  {
707  Int_t relId[4];
708  fPHOSGeo->AbsToRelNumbering(cellList[iCell],relId);
709 
710  irow = relId[2];
711  icol = relId[3];
712  imod = relId[0]-1;
713 
714  if ( fPHOSBadChannelMap->GetEntries() <= imod ) continue;
715 
716  //printf("PHOS bad channels imod %d, icol %d, irow %d\n",imod, irow, icol);
717  if ( GetPHOSChannelStatus(imod, irow, icol) ) return kTRUE;
718  }
719  else return kFALSE;
720  } // cell cluster loop
721 
722  return kFALSE;
723 }
724 
725 //_______________________________________________________________
727 //_______________________________________________________________
729 {
730  clus->SetE(fEMCALRecoUtils->CorrectClusterEnergyLinearity(clus));
731 }
732 
733 //________________________________________________________________________________________
735 //________________________________________________________________________________________
736 Int_t AliCalorimeterUtils::GetMaxEnergyCell(AliVCaloCells* cells, const AliVCluster* clu,
737  Float_t & clusterFraction) const
738 {
739  if( !clu || !cells )
740  {
741  AliInfo("Cluster or cells pointer is null!");
742  return -1;
743  }
744 
745  Double_t eMax =-1.;
746  Double_t eTot = 0.;
747  Double_t eCell =-1.;
748  Float_t fraction = 1.;
749  Float_t recalFactor = 1.;
750  Int_t cellAbsId =-1 , absId =-1 ;
751  Int_t iSupMod =-1 , ieta =-1 , iphi = -1, iRCU = -1;
752 
753  Int_t calo = kEMCAL;
754  if(clu->IsPHOS()) calo = kPHOS ;
755 
756  for (Int_t iDig=0; iDig< clu->GetNCells(); iDig++)
757  {
758  cellAbsId = clu->GetCellAbsId(iDig);
759 
760  fraction = clu->GetCellAmplitudeFraction(iDig);
761  if(fraction < 1e-4) fraction = 1.; // in case unfolding is off
762 
763  iSupMod = GetModuleNumberCellIndexes(cellAbsId, calo, ieta, iphi, iRCU);
764 
765  if(IsRecalibrationOn())
766  {
767  if(calo == kEMCAL) recalFactor = GetEMCALChannelRecalibrationFactor(iSupMod,ieta,iphi);
768  else recalFactor = GetPHOSChannelRecalibrationFactor (iSupMod,iphi,ieta);
769  }
770 
771  eCell = cells->GetCellAmplitude(cellAbsId)*fraction*recalFactor;
772 
773  if(eCell > eMax)
774  {
775  eMax = eCell;
776  absId = cellAbsId;
777  }
778 
779  eTot+=eCell;
780 
781  }// cell loop
782 
783  if(eTot > 0.1)
784  clusterFraction = (eTot-eMax)/eTot; //Do not use cluster energy in case it was corrected for non linearity.
785  else
786  clusterFraction =1.;
787 
788  return absId;
789 }
790 
791 //___________________________________________________________________________________
795 //___________________________________________________________________________________
796 AliVTrack * AliCalorimeterUtils::GetMatchedTrack(AliVCluster* cluster,
797  AliVEvent* event, Int_t index) const
798 {
799  AliVTrack *track = 0x0;
800 
801  //
802  // EMCAL case only when matching is recalculated
803  //
804  if(cluster->IsEMCAL() && IsRecalculationOfClusterTrackMatchingOn())
805  {
806  Int_t trackIndex = fEMCALRecoUtils->GetMatchedTrackIndex(cluster->GetID());
807  //printf("track index %d, cluster ID %d \n ",trackIndex,cluster->GetID());
808 
809  if(trackIndex < 0 )
810  AliInfo(Form("Wrong track index %d, from recalculation", trackIndex));
811  else
812  track = dynamic_cast<AliVTrack*> (event->GetTrack(trackIndex));
813 
814  return track ;
815  }
816 
817  //
818  // Normal case, get info from ESD or AOD
819  //
820 
821  // No tracks matched
822  if( cluster->GetNTracksMatched() < 1 ) return 0x0;
823 
824  // At least one match
825  Int_t iTrack = 0; // only one match for AODs with index 0.
826 
827  // ESDs
828  if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
829  {
830  if( index >= 0 ) iTrack = index;
831  else iTrack = cluster->GetTrackMatchedIndex();
832 
833  track = dynamic_cast<AliVTrack*> ( event->GetTrack(iTrack) );
834  }
835  else // AODs
836  {
837  if( index > 0 ) iTrack = index;
838 
839  track = dynamic_cast<AliVTrack*>( cluster->GetTrackMatched(iTrack) );
840  }
841 
842  return track ;
843 }
844 
849 Float_t AliCalorimeterUtils::GetMCECellClusFracCorrection(Float_t eCell, Float_t eCluster) const
850 {
851  if( eCluster <= 0 || eCluster < eCell )
852  {
853  AliWarning(Form("Bad values eCell=%f, eCluster %f",eCell,eCluster));
854  return 1;
855  }
856 
857  Float_t frac = eCell / eCluster;
858 
859  Float_t correction = fMCECellClusFracCorrParam[0] +
860  TMath::Exp( frac*fMCECellClusFracCorrParam[2]+fMCECellClusFracCorrParam[1] ) +
861  fMCECellClusFracCorrParam[3]/TMath::Sqrt(frac);
862 
863 // printf("AliCalorimeterUtils::GetMCECellClusFracCorrection(eCell=%f, eCluster %f, frac %f) = %f\n",eCell, eCluster, frac, correction);
864 // printf("\t %2.2f + TMath::Exp( %2.3f*%2.2f + %2.2f ) + %2.2f/TMath::Sqrt(%2.3f)) = %f\n",
865 // fMCECellClusFracCorrParam[0],frac,fMCECellClusFracCorrParam[2],fMCECellClusFracCorrParam[1],fMCECellClusFracCorrParam[3], frac, correction);
866 
867  return correction;
868 }
869 
870 //_____________________________________________________________________________________________________
872 //_____________________________________________________________________________________________________
873 Int_t AliCalorimeterUtils::GetModuleNumber(AliAODPWG4Particle * particle, AliVEvent * inputEvent) const
874 {
875  Int_t absId = -1;
876 
877  if(particle->GetDetectorTag()==kEMCAL)
878  {
879  fEMCALGeo->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(), absId);
880 
881  AliDebug(2,Form("EMCAL: cluster eta %f, phi %f, absid %d, SuperModule %d",
882  particle->Eta(), particle->Phi()*TMath::RadToDeg(),absId, fEMCALGeo->GetSuperModuleNumber(absId)));
883 
884  return fEMCALGeo->GetSuperModuleNumber(absId) ;
885  } // EMCAL
886  else if ( particle->GetDetectorTag() == kPHOS )
887  {
888  // In case we use the MC reader, the input are TParticles,
889  // in this case use the corresponing method in PHOS Geometry to get the particle.
890  if(strcmp(inputEvent->ClassName(), "AliMCEvent") == 0 )
891  {
892  Int_t mod =-1;
893  Double_t z = 0., x=0.;
894  TParticle* primary = 0x0;
895  AliStack * stack = ((AliMCEvent*)inputEvent)->Stack();
896 
897  if(stack)
898  {
899  primary = stack->Particle(particle->GetCaloLabel(0));
900  }
901  else
902  {
903  AliFatal("Stack not available, stop!");
904  }
905 
906  if(primary)
907  {
908  fPHOSGeo->ImpactOnEmc(primary,mod,z,x) ;
909  }
910  else
911  {
912  AliFatal("Primary not available, stop!");
913  }
914  return mod;
915  }
916  // Input are ESDs or AODs, get the PHOS module number like this.
917  else
918  {
919  //FIXME
920  //AliVCluster *cluster = inputEvent->GetCaloCluster(particle->GetCaloLabel(0));
921  //return GetModuleNumber(cluster);
922  //MEFIX
923  return -1;
924  }
925  } // PHOS
926 
927  return -1;
928 }
929 
930 //_____________________________________________________________________
932 //_____________________________________________________________________
933 Int_t AliCalorimeterUtils::GetModuleNumber(AliVCluster * cluster) const
934 {
935  if(!cluster)
936  {
937  AliDebug(1,"AliCalorimeterUtils::GetModuleNumber() - NUL Cluster, please check!!!");
938 
939  return -1;
940  }
941 
942  if ( cluster->GetNCells() <= 0 ) return -1;
943 
944  Int_t absId = cluster->GetCellAbsId(0);
945 
946  if ( absId < 0 ) return -1;
947 
948  if( cluster->IsEMCAL() )
949  {
950  AliDebug(2,Form("EMCAL absid %d, SuperModule %d",absId, fEMCALGeo->GetSuperModuleNumber(absId)));
951 
952  return fEMCALGeo->GetSuperModuleNumber(absId) ;
953  } // EMCAL
954  else if ( cluster->IsPHOS() )
955  {
956  Int_t relId[4];
957  fPHOSGeo->AbsToRelNumbering(absId,relId);
958 
959  AliDebug(2,Form("PHOS absid %d Module %d",absId, relId[0]-1));
960 
961  return relId[0]-1;
962  } // PHOS
963 
964  return -1;
965 }
966 
967 //___________________________________________________________________________________________________
969 //___________________________________________________________________________________________________
971  Int_t & icol, Int_t & irow, Int_t & iRCU) const
972 {
973  Int_t imod = -1;
974 
975  if ( absId < 0 ) return -1 ;
976 
977  if ( calo == kEMCAL )
978  {
979  Int_t iTower = -1, iIphi = -1, iIeta = -1;
980  fEMCALGeo->GetCellIndex(absId,imod,iTower,iIphi,iIeta);
981  fEMCALGeo->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,irow,icol);
982 
983  if(imod < 0 || irow < 0 || icol < 0 )
984  {
985  AliFatal(Form("Negative value for super module: %d, or cell icol: %d, or cell irow: %d, check EMCAL geometry name",imod,icol,irow));
986  }
987 
988  // In case of DCal C side, shift columns to match offline/online numbering
989  // when calculating the DDL for Run2
990  // See AliEMCALRawUtils::Digits2Raw and Raw2Digits.
991  Int_t ico2 = icol ;
992  if ( imod == 13 || imod == 15 || imod == 17 ) ico2 += 16;
993 
994  // RCU / DDL
995  if(imod < 10 || (imod > 11 && imod < 18)) // (EMCAL Full || DCAL 2/3)
996  {
997  // RCU0 / DDL0
998  if ( 0 <= irow && irow < 8 ) iRCU = 0; // first cable row
999  else if ( 8 <= irow && irow < 16 &&
1000  0 <= ico2 && ico2 < 24 ) iRCU = 0; // first half;
1001  //second cable row
1002 
1003  // RCU1 / DDL1
1004  else if ( 8 <= irow && irow < 16 &&
1005  24 <= ico2 && ico2 < 48 ) iRCU = 1; // second half;
1006  //second cable row
1007  else if ( 16 <= irow && irow < 24 ) iRCU = 1; // third cable row
1008 
1009  if ( imod%2 == 1 ) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same
1010  }
1011  else
1012  {
1013  // 1/3 SM have one single SRU, just assign RCU/DDL 0
1014  iRCU = 0 ;
1015  }
1016 
1017  if ( iRCU < 0 )
1018  AliFatal(Form("Wrong EMCAL RCU number = %d", iRCU));
1019 
1020  return imod ;
1021  } // EMCAL
1022  else // PHOS
1023  {
1024  Int_t relId[4];
1025  fPHOSGeo->AbsToRelNumbering(absId,relId);
1026 
1027  irow = relId[2];
1028  icol = relId[3];
1029  imod = relId[0]-1;
1030  iRCU= (Int_t)(relId[2]-1)/16 ;
1031 
1032  //Int_t iBranch= (Int_t)(relid[3]-1)/28 ; //0 to 1
1033 
1034  if ( iRCU >= 4 )
1035  AliFatal(Form("Wrong PHOS RCU number = %d", iRCU));
1036 
1037  return imod;
1038  } // PHOS
1039 
1040  return -1;
1041 }
1042 
1043 //___________________________________________________________________________________________
1045 //___________________________________________________________________________________________
1046 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells)
1047 {
1048  const Int_t nc = cluster->GetNCells();
1049 
1050  Int_t absIdList[nc];
1051  Float_t maxEList[nc];
1052 
1053  Int_t nMax = GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList);
1054 
1055  return nMax;
1056 }
1057 
1058 //___________________________________________________________________________________________
1060 //___________________________________________________________________________________________
1061 Int_t AliCalorimeterUtils::GetNumberOfLocalMaxima(AliVCluster* cluster, AliVCaloCells* cells,
1062  Int_t *absIdList, Float_t *maxEList)
1063 {
1064  Int_t iDigitN = 0 ;
1065  Int_t iDigit = 0 ;
1066  Int_t absId1 = -1 ;
1067  Int_t absId2 = -1 ;
1068  const Int_t nCells = cluster->GetNCells();
1069 
1070  Float_t eCluster = RecalibrateClusterEnergy(cluster, cells);// recalculate cluster energy, avoid non lin correction.
1071 
1072  Float_t simuTotWeight = 0;
1074  {
1075  simuTotWeight = RecalibrateClusterEnergyWeightCell(cluster, cells,eCluster);// same but apply a weight
1076  simuTotWeight/= eCluster;
1077  }
1078 
1079  Int_t calorimeter = kEMCAL;
1080  if(!cluster->IsEMCAL()) calorimeter = kPHOS;
1081 
1082  //printf("cluster : ncells %d \n",nCells);
1083 
1084  Float_t emax = 0;
1085  Int_t idmax =-1;
1086  for(iDigit = 0; iDigit < nCells ; iDigit++)
1087  {
1088  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit] ;
1089  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1090  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1091 
1093  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1094 
1095  if( en > emax )
1096  {
1097  emax = en ;
1098  idmax = absIdList[iDigit] ;
1099  }
1100  //Int_t icol = -1, irow = -1, iRCU = -1;
1101  //Int_t sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1102  //printf("\t cell %d, id %d, sm %d, col %d, row %d, e %f\n", iDigit, absIdList[iDigit], sm, icol, irow, en );
1103  }
1104 
1105  for(iDigit = 0 ; iDigit < nCells; iDigit++)
1106  {
1107  if( absIdList[iDigit] >= 0 )
1108  {
1109  absId1 = cluster->GetCellsAbsId()[iDigit];
1110 
1111  Float_t en1 = cells->GetCellAmplitude(absId1);
1112  RecalibrateCellAmplitude(en1,calorimeter,absId1);
1113 
1115  en1*=GetMCECellClusFracCorrection(en1,eCluster)/simuTotWeight;
1116 
1117  //printf("%d : absIDi %d, E %f\n",iDigit, absId1,en1);
1118 
1119  for(iDigitN = 0; iDigitN < nCells; iDigitN++)
1120  {
1121  absId2 = cluster->GetCellsAbsId()[iDigitN] ;
1122 
1123  if(absId2==-1 || absId2==absId1) continue;
1124 
1125  //printf("\t %d : absIDj %d\n",iDigitN, absId2);
1126 
1127  Float_t en2 = cells->GetCellAmplitude(absId2);
1128  RecalibrateCellAmplitude(en2,calorimeter,absId2);
1129 
1131  en2*=GetMCECellClusFracCorrection(en2,eCluster)/simuTotWeight;
1132 
1133  //printf("\t %d : absIDj %d, E %f\n",iDigitN, absId2,en2);
1134 
1135  if ( AreNeighbours(calorimeter, absId1, absId2) )
1136  {
1137  // printf("\t \t Neighbours \n");
1138  if ( en1 > en2 )
1139  {
1140  absIdList[iDigitN] = -1 ;
1141  //printf("\t \t indexN %d not local max\n",iDigitN);
1142  // but may be digit too is not local max ?
1143  if(en1 < en2 + fLocMaxCutEDiff) {
1144  //printf("\t \t index %d not local max cause locMaxCutEDiff\n",iDigit);
1145  absIdList[iDigit] = -1 ;
1146  }
1147  }
1148  else
1149  {
1150  absIdList[iDigit] = -1 ;
1151  //printf("\t \t index %d not local max\n",iDigitN);
1152  // but may be digitN too is not local max ?
1153  if(en1 > en2 - fLocMaxCutEDiff)
1154  {
1155  absIdList[iDigitN] = -1 ;
1156  //printf("\t \t indexN %d not local max cause locMaxCutEDiff\n",iDigit);
1157  }
1158  }
1159  } // if Are neighbours
1160  //else printf("\t \t NOT Neighbours \n");
1161  } // while digitN
1162  } // slot not empty
1163  } // while digit
1164 
1165  iDigitN = 0 ;
1166  for(iDigit = 0; iDigit < nCells; iDigit++)
1167  {
1168  if( absIdList[iDigit] >= 0 )
1169  {
1170  absIdList[iDigitN] = absIdList[iDigit] ;
1171 
1172  Float_t en = cells->GetCellAmplitude(absIdList[iDigit]);
1173  RecalibrateCellAmplitude(en,calorimeter,absIdList[iDigit]);
1174 
1176  en*=GetMCECellClusFracCorrection(en,eCluster)/simuTotWeight;
1177 
1178  if(en < fLocMaxCutE) continue; // Maxima only with seed energy at least
1179 
1180  maxEList[iDigitN] = en ;
1181 
1182  //printf("Local max %d, id %d, en %f\n", iDigit,absIdList[iDigitN],en);
1183  iDigitN++ ;
1184  }
1185  }
1186 
1187  if ( iDigitN == 0 )
1188  {
1189  AliDebug(1,Form("No local maxima found, assign highest energy cell as maxima, id %d, en cell %2.2f, en cluster %2.2f",
1190  idmax,emax,cluster->E()));
1191  iDigitN = 1 ;
1192  maxEList[0] = emax ;
1193  absIdList[0] = idmax ;
1194  }
1195 
1196 
1197  AliDebug(1,Form("In cluster E %2.2f (wth non lin. %2.2f), M02 %2.2f, M20 %2.2f, N maxima %d",
1198  cluster->E(),eCluster, cluster->GetM02(),cluster->GetM20(), iDigitN));
1199 
1200 // if(fDebug > 1) for(Int_t imax = 0; imax < iDigitN; imax++)
1201 // {
1202 // printf(" \t i %d, absId %d, Ecell %f\n",imax,absIdList[imax],maxEList[imax]);
1203 // }
1204 
1205  return iDigitN ;
1206 }
1207 
1208 //____________________________________
1210 //____________________________________
1212 {
1213  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1214  {
1215  AliError("AliCalorimeterUtils::GetPass() - Pointer to tree = 0, returning null");
1216  return TString("");
1217  }
1218 
1219  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1220  {
1221  AliError("AliCalorimeterUtils::GetPass() - Null pointer input file, returning null");
1222  return TString("");
1223  }
1224 
1225  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1226  if (pass.Contains("ass1")) return TString("pass1");
1227  else if (pass.Contains("ass2")) return TString("pass2");
1228  else if (pass.Contains("ass3")) return TString("pass3");
1229  else if (pass.Contains("ass4")) return TString("pass4");
1230  else if (pass.Contains("ass5")) return TString("pass5");
1231  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1232  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1233  {
1234  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1235  return TString("pass1");
1236  }
1237 
1238  // No condition fullfilled, give a default value
1239  AliInfo("Pass number string not found");
1240  return TString("");
1241 }
1242 
1243 //________________________________________
1245 //________________________________________
1247 {
1248  fEMCALGeoName = "";
1249  fPHOSGeoName = "";
1250 
1251  fEMCALGeoMatrixSet = kFALSE;
1252  fPHOSGeoMatrixSet = kFALSE;
1253 
1254  fRemoveBadChannels = kFALSE;
1255 
1257 
1258  fLocMaxCutE = 0.1 ;
1259  fLocMaxCutEDiff = 0.0 ;
1260 
1261  // fMaskCellColumns = new Int_t[fNMaskCellColumns];
1262  // fMaskCellColumns[0] = 6 ; fMaskCellColumns[1] = 7 ; fMaskCellColumns[2] = 8 ;
1263  // fMaskCellColumns[3] = 35; fMaskCellColumns[4] = 36; fMaskCellColumns[5] = 37;
1264  // fMaskCellColumns[6] = 12+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[7] = 13+AliEMCALGeoParams::fgkEMCALCols;
1265  // fMaskCellColumns[8] = 40+AliEMCALGeoParams::fgkEMCALCols; fMaskCellColumns[9] = 41+AliEMCALGeoParams::fgkEMCALCols;
1266  // fMaskCellColumns[10]= 42+AliEMCALGeoParams::fgkEMCALCols;
1267 
1268  fOADBSet = kFALSE;
1269  fOADBForEMCAL = kTRUE ;
1270  fOADBForPHOS = kFALSE;
1271 
1272  fOADBFilePathEMCAL = "$ALICE_PHYSICS/OADB/EMCAL" ;
1273  fOADBFilePathPHOS = "$ALICE_PHYSICS/OADB/PHOS" ;
1274 
1275  fImportGeometryFromFile = kTRUE;
1277 
1278  fNSuperModulesUsed = 22;
1279 
1280  fMCECellClusFracCorrParam[0] = 0.78;
1281  fMCECellClusFracCorrParam[1] =-1.8;
1282  fMCECellClusFracCorrParam[2] =-6.3;
1283  fMCECellClusFracCorrParam[3] = 0.014;
1284 }
1285 
1286 //_____________________________________________________
1288 //_____________________________________________________
1290 {
1291  AliDebug(1,"Init bad channel map");
1292 
1293  // In order to avoid rewriting the same histograms
1294  Bool_t oldStatus = TH1::AddDirectoryStatus();
1295  TH1::AddDirectory(kFALSE);
1296 
1297  fPHOSBadChannelMap = new TObjArray(5);
1298  for (int i = 0; i < 5; i++)fPHOSBadChannelMap->Add(new TH2I(Form("PHOS_BadMap_mod%d",i),
1299  Form("PHOS_BadMap_mod%d",i),
1300  64, 0, 64, 56, 0, 56));
1301 
1302  fPHOSBadChannelMap->SetOwner(kTRUE);
1303  fPHOSBadChannelMap->Compress();
1304 
1305  //In order to avoid rewriting the same histograms
1306  TH1::AddDirectory(oldStatus);
1307 }
1308 
1309 //______________________________________________________
1311 //______________________________________________________
1313 {
1314  AliDebug(1,"Init recalibration map");
1315 
1316  // In order to avoid rewriting the same histograms
1317  Bool_t oldStatus = TH1::AddDirectoryStatus();
1318  TH1::AddDirectory(kFALSE);
1319 
1320  fPHOSRecalibrationFactors = new TObjArray(5);
1321  for (int i = 0; i < 5; i++)
1322  {
1323  fPHOSRecalibrationFactors->Add(new TH2F(Form("PHOSRecalFactors_Mod%d",i),
1324  Form("PHOSRecalFactors_Mod%d",i),
1325  64, 0, 64, 56, 0, 56));
1326  }
1327 
1328  // Init the histograms with 1
1329  for (Int_t m = 0; m < 5; m++)
1330  {
1331  for (Int_t i = 0; i < 56; i++)
1332  {
1333  for (Int_t j = 0; j < 64; j++)
1334  {
1336  }
1337  }
1338  }
1339 
1340  fPHOSRecalibrationFactors->SetOwner(kTRUE);
1341  fPHOSRecalibrationFactors->Compress();
1342 
1343  // In order to avoid rewriting the same histograms
1344  TH1::AddDirectory(oldStatus);
1345 }
1346 
1351 {
1352  if (fEMCALGeo) return;
1353 
1354  AliDebug(1,Form(" for run=%d",fRunNumber));
1355 
1356  if(fEMCALGeoName=="")
1357  {
1358  fEMCALGeo = AliEMCALGeometry::GetInstanceFromRunNumber(fRunNumber);
1359  AliInfo(Form("Get EMCAL geometry name to <%s> for run %d",fEMCALGeo->GetName(),fRunNumber));
1360  }
1361  else
1362  {
1363  fEMCALGeo = AliEMCALGeometry::GetInstance(fEMCALGeoName);
1364  AliInfo(Form("Set EMCAL geometry name to <%s>",fEMCALGeoName.Data()));
1365  }
1366 
1367  // Init geometry, I do not like much to do it like this ...
1368  if(fImportGeometryFromFile && !gGeoManager)
1369  {
1370  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1371  {
1372  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1373  if (fRunNumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1374  else if(fRunNumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1375  else if(fRunNumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1376  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >= 2015
1377  }
1378 
1379  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1380 
1381  TGeoManager::Import(fImportGeometryFilePath) ; // default need file "geometry.root" in local dir!!!!
1382  }
1383  else if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1384 }
1385 
1390 {
1391  if (fPHOSGeo) return;
1392 
1393  AliDebug(1,Form(" for run=%d",fRunNumber));
1394 
1395  if(fPHOSGeoName=="") fPHOSGeoName = "PHOSgeo";
1396 
1397  fPHOSGeo = new AliPHOSGeoUtils(fPHOSGeoName);
1398 
1399  //if (!gGeoManager) AliInfo("Careful!, gGeoManager not loaded, load misalign matrices");
1400 }
1401 
1402 //______________________________________________________________________________________________
1403 // Check that a MC ESD is in the calorimeter acceptance
1404 //______________________________________________________________________________________________
1405 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(Int_t calo, TParticle* particle)
1406 {
1407  if(!particle || (calo!=kEMCAL && calo!=kPHOS)) return kFALSE ;
1408 
1409  if( (!IsPHOSGeoMatrixSet () && calo == kPHOS ) ||
1410  (!IsEMCALGeoMatrixSet() && calo == kEMCAL) )
1411  {
1412  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1413  return kFALSE ;
1414  }
1415 
1416  if(calo == kPHOS )
1417  {
1418  Int_t mod = 0 ;
1419  Double_t x = 0, z = 0 ;
1420  return GetPHOSGeometry()->ImpactOnEmc( particle, mod, z, x);
1421  }
1422  else if(calo == kEMCAL)
1423  {
1424  Int_t absID = 0 ;
1425  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),particle->Phi(),absID);
1426  if(ok)
1427  {
1428  Int_t icol = -1, irow = -1, iRCU = -1;
1429  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1430  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1431  if(status > 0) ok = kFALSE;
1432  }
1433 
1434  return ok ;
1435  }
1436 
1437  return kFALSE ;
1438 }
1439 
1440 //______________________________________________________________________________________________________
1442 //______________________________________________________________________________________________________
1443 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(Int_t calo, AliAODMCParticle* particle)
1444 {
1445  if(!particle || (calo!=kEMCAL && calo!=kPHOS)) return kFALSE ;
1446 
1447  if( (!IsPHOSGeoMatrixSet () && calo == kPHOS ) ||
1448  (!IsEMCALGeoMatrixSet() && calo == kEMCAL) )
1449  {
1450  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1451  return kFALSE ;
1452  }
1453 
1454  Float_t phi = particle->Phi();
1455  if(phi < 0) phi+=TMath::TwoPi();
1456 
1457  if(calo == kPHOS )
1458  {
1459  Int_t mod = 0 ;
1460  Double_t x = 0, z = 0 ;
1461  Double_t vtx[]={ particle->Xv(), particle->Yv(), particle->Zv() } ;
1462  return GetPHOSGeometry()->ImpactOnEmc(vtx, particle->Theta(), phi, mod, z, x) ;
1463  }
1464  else if(calo == kEMCAL)
1465  {
1466  Int_t absID = 0 ;
1467  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(particle->Eta(),phi,absID);
1468  if(ok)
1469  {
1470  Int_t icol = -1, irow = -1, iRCU = -1;
1471  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1472  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1473  if(status > 0) ok = kFALSE;
1474  }
1475 
1476  return ok ;
1477  }
1478 
1479  return kFALSE ;
1480 }
1481 
1482 //_____________________________________________________________________________________________________
1483 // Check that a TLorentzVector is in the calorimeter acceptance, give the cell number where it hit.
1484 //_____________________________________________________________________________________________________
1485 Bool_t AliCalorimeterUtils::IsMCParticleInCalorimeterAcceptance(Int_t calo, Float_t eta, Float_t theta,
1486  Float_t phiOrg, Int_t & absID)
1487 {
1488  if(calo!=kEMCAL && calo!=kPHOS) return kFALSE ;
1489 
1490  if( (!IsPHOSGeoMatrixSet () && calo == kPHOS ) ||
1491  (!IsEMCALGeoMatrixSet() && calo == kEMCAL) )
1492  {
1493  AliFatal(Form("Careful Geo Matrix for calo <%d> is not set, use AliFidutialCut instead",calo));
1494  return kFALSE ;
1495  }
1496 
1497  Float_t phi = phiOrg;
1498  if(phi < 0) phi+=TMath::TwoPi();
1499 
1500  if(calo == kPHOS )
1501  {
1502  Int_t mod = 0 ;
1503  Double_t x = 0, z = 0 ;
1504  Double_t vtx[]={0,0,0} ;
1505  return GetPHOSGeometry()->ImpactOnEmc(vtx, theta, phi, mod, z, x) ;
1506  }
1507  else if(calo == kEMCAL)
1508  {
1509  Bool_t ok = GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(eta,phi,absID);
1510  if(ok)
1511  {
1512  Int_t icol = -1, irow = -1, iRCU = -1;
1513  Int_t nModule = GetModuleNumberCellIndexes(absID,calo, icol, irow, iRCU);
1514  Int_t status = GetEMCALChannelStatus(nModule,icol,irow);
1515  if(status > 0) ok = kFALSE;
1516  }
1517 
1518  return ok ;
1519  }
1520 
1521  return kFALSE ;
1522 }
1523 
1524 //_______________________________________________________________________
1527 //_______________________________________________________________________
1528 Bool_t AliCalorimeterUtils::MaskFrameCluster(Int_t iSM, Int_t ieta) const
1529 {
1530  Int_t icol = ieta;
1531  if ( iSM%2 ) icol+=48; // Impair SM, shift index [0-47] to [48-96]
1532 
1534  {
1535  for (Int_t imask = 0; imask < fNMaskCellColumns; imask++)
1536  {
1537  if(icol==fMaskCellColumns[imask]) return kTRUE;
1538  }
1539  }
1540 
1541  return kFALSE;
1542 }
1543 
1544 //_________________________________________________________
1546 //_________________________________________________________
1547 void AliCalorimeterUtils::Print(const Option_t * opt) const
1548 {
1549  if(! opt)
1550  return;
1551 
1552  printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
1553  printf("Remove Clusters with bad channels? %d\n",fRemoveBadChannels);
1554  printf("Remove Clusters with max cell at less than %d cells from EMCAL border and %d cells from PHOS border\n",
1555  fEMCALRecoUtils->GetNumberOfCellsFromEMCALBorder(), fNCellsFromPHOSBorder);
1556  if(fEMCALRecoUtils->IsEMCALNoBorderAtEta0()) printf("Do not remove EMCAL clusters at Eta = 0\n");
1557  printf("Recalibrate Clusters? %d, run by run %d\n",fRecalibration,fRunDependentCorrection);
1558  printf("Recalculate Clusters Position? %d\n",fRecalculatePosition);
1559  printf("Recalculate Clusters Energy? %d\n",fCorrectELinearity);
1560  printf("Matching criteria: dR < %2.2f[cm], dZ < %2.2f[cm]\n",fCutR,fCutZ);
1561 
1562  printf("Loc. Max. E > %2.2f\n", fLocMaxCutE);
1563  printf("Loc. Max. E Diff > %2.2f\n", fLocMaxCutEDiff);
1564 
1565  printf(" \n") ;
1566 }
1567 
1568 //_____________________________________________________________________________________________
1570 //_____________________________________________________________________________________________
1571 void AliCalorimeterUtils::RecalibrateCellAmplitude(Float_t & amp, Int_t calo, Int_t id) const
1572 {
1573  Int_t icol = -1; Int_t irow = -1; Int_t iRCU = -1;
1574  Int_t nModule = GetModuleNumberCellIndexes(id,calo, icol, irow, iRCU);
1575 
1576  if (IsRecalibrationOn())
1577  {
1578  if(calo == kPHOS)
1579  {
1580  amp *= GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
1581  }
1582  else
1583  {
1584  amp *= GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
1585  }
1586  }
1587 }
1588 
1589 //____________________________________________________________________________________________________
1591 //____________________________________________________________________________________________________
1592 void AliCalorimeterUtils::RecalibrateCellTime(Double_t & time, Int_t calo, Int_t id, Int_t bc) const
1593 {
1594  if ( calo == kEMCAL && GetEMCALRecoUtils()->IsTimeRecalibrationOn() )
1595  {
1596  GetEMCALRecoUtils()->RecalibrateCellTime(id,bc,time);
1597  }
1598 }
1599 
1600 //__________________________________________________________________________
1602 //__________________________________________________________________________
1603 Float_t AliCalorimeterUtils::RecalibrateClusterEnergy(AliVCluster * cluster,
1604  AliVCaloCells * cells)
1605 {
1606  // Initialize some used variables
1607  Float_t frac = 0., energy = 0.;
1608 
1609  if(cells)
1610  {
1611  //Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1612 
1613  UShort_t * index = cluster->GetCellsAbsId() ;
1614  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1615 
1616  Int_t ncells = cluster->GetNCells();
1617 
1618  Int_t calo = kEMCAL;
1619  if(cluster->IsPHOS()) calo = kPHOS ;
1620 
1621  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1622  for(Int_t icell = 0; icell < ncells; icell++)
1623  {
1624  Int_t absId = index[icell];
1625 
1626  frac = fraction[icell];
1627  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1628 
1629  Float_t amp = cells->GetCellAmplitude(absId);
1630  RecalibrateCellAmplitude(amp,calo, absId);
1631 
1632  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy: before cal %f; after cal %f",
1633  calo,frac,cells->GetCellAmplitude(absId),amp));
1634 
1635  energy += amp*frac;
1636  }
1637 
1638  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
1639 
1640  } // cells available
1641  else
1642  {
1643  AliFatal("Cells pointer does not exist!");
1644  }
1645 
1646  return energy;
1647 }
1648 
1649 //_______________________________________________________________________________________________________
1652 //_______________________________________________________________________________________________________
1654  AliVCaloCells * cells, Float_t energyOrg)
1655 {
1656  //Initialize some used variables
1657  Float_t frac = 0., energy = 0.;
1658 
1659  if(cells)
1660  {
1661  // Get the cluster number of cells and list of absId, check what kind of cluster do we have.
1662 
1663  UShort_t * index = cluster->GetCellsAbsId() ;
1664  Double_t * fraction = cluster->GetCellsAmplitudeFraction() ;
1665 
1666  Int_t ncells = cluster->GetNCells();
1667 
1668  Int_t calo = kEMCAL;
1669  if(cluster->IsPHOS()) calo = kPHOS ;
1670 
1671  // Loop on the cells, get the cell amplitude and recalibration factor, multiply and and to the new energy
1672  for(Int_t icell = 0; icell < ncells; icell++)
1673  {
1674  Int_t absId = index[icell];
1675 
1676  frac = fraction[icell];
1677  if(frac < 1e-3) frac = 1; //in case of EMCAL, this is set as 0, not used.
1678 
1679  Float_t amp = cells->GetCellAmplitude(absId);
1680  RecalibrateCellAmplitude(amp,calo, absId);
1681 
1682  amp*=GetMCECellClusFracCorrection(amp,energyOrg);
1683 
1684  AliDebug(2,Form("Recalibrate cell: calo <%d>, cell fraction %f, cell energy %f",
1685  calo,frac,cells->GetCellAmplitude(absId)));
1686 
1687  energy += amp*frac;
1688  }
1689 
1690  AliDebug(1,Form("Energy before %f, after %f",cluster->E(),energy));
1691  } // cells available
1692  else
1693  {
1694  AliFatal("Cells pointer does not exist!");
1695  }
1696 
1697  return energy;
1698 }
1699 
1700 //__________________________________________________________________________________________
1702 //__________________________________________________________________________________________
1703 void AliCalorimeterUtils::RecalculateClusterPosition(AliVCaloCells* cells, AliVCluster* clu)
1704 {
1705  fEMCALRecoUtils->RecalculateClusterPosition((AliEMCALGeometry*)fEMCALGeo, cells,clu);
1706 }
1707 
1708 //________________________________________________________________________________
1710 //________________________________________________________________________________
1712  TObjArray* clusterArray)
1713 {
1714  if (fRecalculateMatching)
1715  {
1716  fEMCALRecoUtils->FindMatches(event,clusterArray,fEMCALGeo) ;
1717  //AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
1718  //if(esdevent){
1719  // fEMCALRecoUtils->SetClusterMatchedToTrack(esdevent) ;
1720  // fEMCALRecoUtils->SetTracksMatchedToCluster(esdevent) ;
1721  //}
1722  }
1723 }
1724 
1725 //___________________________________________________________________________
1738 //___________________________________________________________________________
1739 void AliCalorimeterUtils::SplitEnergy(Int_t absId1, Int_t absId2,
1740  AliVCluster* cluster,
1741  AliVCaloCells* cells,
1742  //Float_t & e1, Float_t & e2,
1743  AliAODCaloCluster* cluster1,
1744  AliAODCaloCluster* cluster2,
1745  Int_t nMax, Int_t eventNumber)
1746 {
1747  TH2F* hClusterMap = 0 ;
1748  TH2F* hClusterLocMax = 0 ;
1749  TH2F* hCluster1 = 0 ;
1750  TH2F* hCluster2 = 0 ;
1751 
1752  if(fPlotCluster)
1753  {
1754  hClusterMap = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
1755  hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
1756  hCluster1 = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
1757  hCluster2 = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
1758  hClusterMap ->SetXTitle("column");
1759  hClusterMap ->SetYTitle("row");
1760  hClusterLocMax ->SetXTitle("column");
1761  hClusterLocMax ->SetYTitle("row");
1762  hCluster1 ->SetXTitle("column");
1763  hCluster1 ->SetYTitle("row");
1764  hCluster2 ->SetXTitle("column");
1765  hCluster2 ->SetYTitle("row");
1766  }
1767 
1768  Int_t calorimeter = kEMCAL;
1769  if(cluster->IsPHOS())
1770  {
1771  calorimeter=kPHOS;
1772  AliWarning("Not supported for PHOS yet");
1773  return;
1774  }
1775 
1776  const Int_t ncells = cluster->GetNCells();
1777  Int_t absIdList[ncells];
1778 
1779  Float_t e1 = 0, e2 = 0 ;
1780  Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;
1781  Float_t eCluster = 0;
1782  Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1;
1783  for(Int_t iDigit = 0; iDigit < ncells; iDigit++ )
1784  {
1785  absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
1786 
1787  Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
1788  RecalibrateCellAmplitude(ec,calorimeter, absIdList[iDigit]);
1789  eCluster+=ec;
1790 
1791  if(fPlotCluster)
1792  {
1793  //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
1794  sm = GetModuleNumberCellIndexes(absIdList[iDigit], calorimeter, icol, irow, iRCU) ;
1795  if(sm > -1 && sm < 12) // just to avoid compilation warning
1796  {
1797  if(icol > maxCol) maxCol = icol;
1798  if(icol < minCol) minCol = icol;
1799  if(irow > maxRow) maxRow = irow;
1800  if(irow < minRow) minRow = irow;
1801  hClusterMap->Fill(icol,irow,ec);
1802  }
1803  }
1804  }
1805 
1806  // Init counters and variables
1807  Int_t ncells1 = 1 ;
1808  UShort_t absIdList1[9] ;
1809  Double_t fracList1 [9] ;
1810  absIdList1[0] = absId1 ;
1811  fracList1 [0] = 1. ;
1812 
1813  Float_t ecell1 = cells->GetCellAmplitude(absId1);
1814  RecalibrateCellAmplitude(ecell1, calorimeter, absId1);
1815  e1 = ecell1;
1816 
1817  Int_t ncells2 = 1 ;
1818  UShort_t absIdList2[9] ;
1819  Double_t fracList2 [9] ;
1820  absIdList2[0] = absId2 ;
1821  fracList2 [0] = 1. ;
1822 
1823  Float_t ecell2 = cells->GetCellAmplitude(absId2);
1824  RecalibrateCellAmplitude(ecell2, calorimeter, absId2);
1825  e2 = ecell2;
1826 
1827  if(fPlotCluster)
1828  {
1829  Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
1830  sm = GetModuleNumberCellIndexes(absId1, calorimeter, icol1, irow1, iRCU) ;
1831  hClusterLocMax->Fill(icol1,irow1,ecell1);
1832  sm = GetModuleNumberCellIndexes(absId2, calorimeter, icol2, irow2, iRCU) ;
1833  hClusterLocMax->Fill(icol2,irow2,ecell2);
1834  }
1835 
1836  // Very rough way to share the cluster energy
1837  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
1838  Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
1839  Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
1840 
1841  for(Int_t iDigit = 0; iDigit < ncells; iDigit++)
1842  {
1843  Int_t absId = absIdList[iDigit];
1844 
1845  if ( absId==absId1 || absId==absId2 || absId < 0 ) continue;
1846 
1847  Float_t ecell = cells->GetCellAmplitude(absId);
1848  RecalibrateCellAmplitude(ecell, calorimeter, absId);
1849 
1850  if(AreNeighbours(calorimeter, absId1,absId ))
1851  {
1852  absIdList1[ncells1]= absId;
1853 
1854  if(AreNeighbours(calorimeter, absId2,absId ))
1855  {
1856  fracList1[ncells1] = shareFraction1;
1857  e1 += ecell*shareFraction1;
1858  }
1859  else
1860  {
1861  fracList1[ncells1] = 1.;
1862  e1 += ecell;
1863  }
1864 
1865  ncells1++;
1866 
1867  } // neigbour to cell1
1868 
1869  if(AreNeighbours(calorimeter, absId2,absId ))
1870  {
1871  absIdList2[ncells2]= absId;
1872 
1873  if(AreNeighbours(calorimeter, absId1,absId ))
1874  {
1875  fracList2[ncells2] = shareFraction2;
1876  e2 += ecell*shareFraction2;
1877  }
1878  else
1879  {
1880  fracList2[ncells2] = 1.;
1881  e2 += ecell;
1882  }
1883 
1884  ncells2++;
1885 
1886  } // neigbour to cell2
1887  }
1888 
1889  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",
1890  nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2));
1891 
1892  cluster1->SetE(e1);
1893  cluster2->SetE(e2);
1894 
1895  cluster1->SetNCells(ncells1);
1896  cluster2->SetNCells(ncells2);
1897 
1898  cluster1->SetCellsAbsId(absIdList1);
1899  cluster2->SetCellsAbsId(absIdList2);
1900 
1901  cluster1->SetCellsAmplitudeFraction(fracList1);
1902  cluster2->SetCellsAmplitudeFraction(fracList2);
1903 
1904  // Correct linearity
1905  CorrectClusterEnergy(cluster1) ;
1906  CorrectClusterEnergy(cluster2) ;
1907 
1908  if(calorimeter==kEMCAL)
1909  {
1910  GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
1911  GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
1912  }
1913 
1914  if(fPlotCluster)
1915  {
1916  //printf("Cells of cluster1: ");
1917  for(Int_t iDigit = 0; iDigit < ncells1; iDigit++ )
1918  {
1919  //printf(" %d ",absIdList1[iDigit]);
1920 
1921  sm = GetModuleNumberCellIndexes(absIdList1[iDigit], calorimeter, icol, irow, iRCU) ;
1922 
1923  Float_t ecell = cells->GetCellAmplitude(absIdList1[iDigit]);
1924  RecalibrateCellAmplitude(ecell, calorimeter, absIdList1[iDigit]);
1925 
1926  if( AreNeighbours(calorimeter, absId2,absIdList1[iDigit]) && absId1!=absIdList1[iDigit])
1927  hCluster1->Fill(icol,irow,ecell*shareFraction1);
1928  else
1929  hCluster1->Fill(icol,irow,ecell);
1930  }
1931 
1932  //printf(" \n ");
1933  //printf("Cells of cluster2: ");
1934 
1935  for(Int_t iDigit = 0; iDigit < ncells2; iDigit++ )
1936  {
1937  //printf(" %d ",absIdList2[iDigit]);
1938 
1939  sm = GetModuleNumberCellIndexes(absIdList2[iDigit], calorimeter, icol, irow, iRCU) ;
1940 
1941  Float_t ecell = cells->GetCellAmplitude(absIdList2[iDigit]);
1942  RecalibrateCellAmplitude(ecell, calorimeter, absIdList2[iDigit]);
1943 
1944  if( AreNeighbours(calorimeter, absId1,absIdList2[iDigit]) && absId2!=absIdList2[iDigit])
1945  hCluster2->Fill(icol,irow,ecell*shareFraction2);
1946  else
1947  hCluster2->Fill(icol,irow,ecell);
1948  }
1949  //printf(" \n ");
1950 
1951  gStyle->SetPadRightMargin(0.1);
1952  gStyle->SetPadLeftMargin(0.1);
1953  gStyle->SetOptStat(0);
1954  gStyle->SetOptFit(000000);
1955 
1956  if(maxCol-minCol > maxRow-minRow)
1957  {
1958  maxRow+= maxCol-minCol;
1959  }
1960  else
1961  {
1962  maxCol+= maxRow-minRow;
1963  }
1964 
1965  TCanvas * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
1966  c->Divide(2,2);
1967  c->cd(1);
1968  gPad->SetGridy();
1969  gPad->SetGridx();
1970  gPad->SetLogz();
1971  hClusterMap ->SetAxisRange(minCol, maxCol,"X");
1972  hClusterMap ->SetAxisRange(minRow, maxRow,"Y");
1973  hClusterMap ->Draw("colz TEXT");
1974  c->cd(2);
1975  gPad->SetGridy();
1976  gPad->SetGridx();
1977  gPad->SetLogz();
1978  hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
1979  hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
1980  hClusterLocMax ->Draw("colz TEXT");
1981  c->cd(3);
1982  gPad->SetGridy();
1983  gPad->SetGridx();
1984  gPad->SetLogz();
1985  hCluster1 ->SetAxisRange(minCol, maxCol,"X");
1986  hCluster1 ->SetAxisRange(minRow, maxRow,"Y");
1987  hCluster1 ->Draw("colz TEXT");
1988  c->cd(4);
1989  gPad->SetGridy();
1990  gPad->SetGridx();
1991  gPad->SetLogz();
1992  hCluster2 ->SetAxisRange(minCol, maxCol,"X");
1993  hCluster2 ->SetAxisRange(minRow, maxRow,"Y");
1994  hCluster2 ->Draw("colz TEXT");
1995 
1996  if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",
1997  eventNumber,cluster->E(),nMax,ncells1,ncells2));
1998 
1999  delete c;
2000  delete hClusterMap;
2001  delete hClusterLocMax;
2002  delete hCluster1;
2003  delete hCluster2;
2004  }
2005 }
2006 
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)
Recalculate track matching.
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
Int_t GetMaxEnergyCell(AliVCaloCells *cells, const AliVCluster *clu, Float_t &fraction) const
For a given CaloCluster, it gets the absId of the cell with maximum energy deposit.
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)
Recalculate EMCAL cluster position.
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.
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...
Bool_t CheckCellFiducialRegion(AliVCluster *cluster, AliVCaloCells *cells) const