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