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