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