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