Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : // --- ROOT system ---
17 : #include <TClonesArray.h>
18 : #include "TGeoManager.h"
19 : #include "TGeoMatrix.h"
20 :
21 : // --- Standard library ---
22 :
23 : // --- AliRoot header files ---
24 : #include "AliEMCALReconstructor.h"
25 :
26 : #include "AliCodeTimer.h"
27 : #include "AliCaloCalibPedestal.h"
28 : #include "AliEMCALCalibData.h"
29 : #include "AliEMCALCalibTime.h"
30 : #include "AliESDEvent.h"
31 : #include "AliESDCaloCluster.h"
32 : #include "AliESDCaloCells.h"
33 : #include "AliESDtrack.h"
34 : #include "AliEMCALLoader.h"
35 : #include "AliEMCALRawUtils.h"
36 : #include "AliEMCALDigit.h"
37 : #include "AliEMCALClusterizerv1.h"
38 : #include "AliEMCALClusterizerv2.h"
39 : #include "AliEMCALClusterizerNxN.h"
40 : #include "AliEMCALRecPoint.h"
41 : #include "AliEMCALPID.h"
42 : #include "AliEMCALRecoUtils.h"
43 : #include "AliRawReader.h"
44 : #include "AliCDBEntry.h"
45 : #include "AliCDBManager.h"
46 : #include "AliEMCALGeometry.h"
47 : #include "AliEMCAL.h"
48 : #include "AliESDVZERO.h"
49 : #include "AliCDBManager.h"
50 : #include "AliRunLoader.h"
51 : #include "AliRun.h"
52 : #include "AliEMCALTriggerData.h"
53 : #include "AliEMCALTriggerElectronics.h"
54 : #include "AliEMCALTriggerDCSConfigDB.h"
55 : #include "AliEMCALTriggerDCSConfig.h"
56 : #include "AliEMCALTriggerData.h"
57 : #include "AliEMCALTriggerRawDigit.h"
58 : #include "AliEMCALTriggerPatch.h"
59 : #include "AliEMCALTriggerTypes.h"
60 :
61 : /// \cond CLASSIMP
62 42 : ClassImp(AliEMCALReconstructor) ;
63 : /// \endcond
64 :
65 : const AliEMCALRecParam* AliEMCALReconstructor::fgkRecParam = 0; // EMCAL rec. parameters
66 : AliEMCALRawUtils* AliEMCALReconstructor::fgRawUtils = 0; // EMCAL raw utilities class
67 : AliEMCALClusterizer* AliEMCALReconstructor::fgClusterizer = 0; // EMCAL clusterizer class
68 : TClonesArray* AliEMCALReconstructor::fgDigitsArr = 0; // list of digits, to be used multiple times
69 : TObjArray* AliEMCALReconstructor::fgClustersArr = 0; // list of clusters, to be used multiple times
70 : TClonesArray* AliEMCALReconstructor::fgTriggerDigits = 0; // list of trigger digits, to be used multiple times
71 : AliEMCALTriggerElectronics *AliEMCALReconstructor::fgTriggerProcessor = 0x0;
72 : TClonesArray *AliEMCALReconstructor::fgTriggerData = 0x0;
73 :
74 : ///
75 : /// Constructor
76 : /// Init all arrays and stuff.
77 : ///
78 : //____________________________________________________________________________
79 2 : AliEMCALReconstructor::AliEMCALReconstructor()
80 2 : : fGeom(0),fCalibData(0),fCalibTime(0),fPedestalData(0), fMatches(0x0)
81 10 : {
82 6 : fgRawUtils = new AliEMCALRawUtils;
83 :
84 2 : AliCDBManager* man = AliCDBManager::Instance();
85 :
86 : //----------------------------------
87 : // Get the geometry, 3 posibilities:
88 : // * To make sure we match with the geometry in a simulation file,
89 : // let's try to get it first.
90 : // * If not, check the run number assigned for this chunk and set the
91 : // geometry depending on the run number
92 : // * If not, take the default geometry
93 :
94 2 : AliRunLoader *rl = AliRunLoader::Instance();
95 4 : if (rl->GetAliRun())
96 : {
97 5 : AliEMCAL * emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
98 3 : if(emcal) fGeom = emcal->GetGeometry();
99 1 : }
100 :
101 2 : if(!fGeom)
102 : {
103 1 : Int_t runNumber = man->GetRun();
104 3 : fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(runNumber);
105 1 : }
106 :
107 2 : if(!fGeom)
108 : {
109 0 : AliWarning(Form("Using default geometry in reconstruction!!!"));
110 0 : fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
111 0 : }
112 :
113 2 : if ( !fGeom ) AliFatal(Form("Could not get geometry!"));
114 8 : else AliInfo (Form("Geometry name: <<%s>>",fGeom->GetName()));
115 :
116 : //---------------------------
117 : // Get energy calibration parameters
118 2 : if(!fCalibData)
119 : {
120 6 : AliCDBEntry *entry = (AliCDBEntry*) man->Get("EMCAL/Calib/Data");
121 4 : if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
122 2 : }
123 :
124 2 : if(!fCalibData)
125 0 : AliFatal("Energy Calibration parameters not found in CDB!");
126 :
127 :
128 : //---------------------------
129 : // Get time calibration parameters if requested
130 2 : if(!fCalibTime)
131 : {
132 6 : AliCDBEntry *entry = (AliCDBEntry*) man->Get("EMCAL/Calib/Time");
133 4 : if (entry) fCalibTime = (AliEMCALCalibTime*) entry->GetObject();
134 2 : }
135 :
136 2 : if(!fCalibTime)
137 0 : AliFatal("Time Calibration parameters not found in CDB!");
138 :
139 : //------------------
140 : // Get bad channels
141 2 : if(!fPedestalData)
142 : {
143 6 : AliCDBEntry *entry = (AliCDBEntry*) man->Get("EMCAL/Calib/Pedestals");
144 4 : if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
145 2 : }
146 :
147 2 : if(!fPedestalData)
148 0 : AliFatal("Dead map not found in CDB!");
149 :
150 : //----------------------------------------------------
151 : // Get trigger parameters and init other trigger stuff
152 2 : AliEMCALTriggerDCSConfigDB* dcsConfigDB = AliEMCALTriggerDCSConfigDB::Instance();
153 :
154 2 : const AliEMCALTriggerDCSConfig* dcsConfig = dcsConfigDB->GetTriggerDCSConfig();
155 :
156 2 : if (!dcsConfig)
157 0 : AliFatal("No Trigger DCS Configuration from OCDB!");
158 :
159 6 : fgTriggerProcessor = new AliEMCALTriggerElectronics( dcsConfig );
160 :
161 4 : int dsize = (fGeom->GetTriggerMappingVersion() == 2) ? 2 : 1;
162 6 : fgTriggerData = new TClonesArray("AliEMCALTriggerData",dsize);
163 10 : for (int i=0;i<dsize;i++) {
164 9 : new((*fgTriggerData)[i]) AliEMCALTriggerData();
165 : }
166 :
167 : //-----------------------------
168 : // Init temporary list of digits
169 6 : fgDigitsArr = new TClonesArray("AliEMCALDigit",1000);
170 6 : fgClustersArr = new TObjArray(1000);
171 :
172 2 : const int kNTRU = fGeom->GetNTotalTRU();
173 6 : fgTriggerDigits = new TClonesArray("AliEMCALTriggerRawDigit", kNTRU * 96);
174 :
175 : //--------------------------
176 : // Init Track matching array
177 6 : fMatches = new TList();
178 2 : fMatches->SetOwner(kTRUE);
179 4 : }
180 :
181 : ///
182 : /// Destructor
183 : ///
184 : //____________________________________________________________________________
185 : AliEMCALReconstructor::~AliEMCALReconstructor()
186 12 : {
187 : ////RS if(fGeom) delete fGeom;
188 :
189 : //No need to delete, recovered from OCDB
190 : //if(fCalibData) delete fCalibData;
191 : //if(fCalibTime) delete fCalibTime;
192 : //if(fPedestalData) delete fPedestalData;
193 :
194 4 : if(fgDigitsArr) fgDigitsArr->Clear("C");
195 4 : delete fgDigitsArr;
196 2 : fgDigitsArr = 0;
197 :
198 4 : if(fgClustersArr) fgClustersArr->Clear();
199 4 : delete fgClustersArr;
200 2 : fgClustersArr = 0;
201 :
202 4 : if(fgTriggerDigits) fgTriggerDigits->Clear();
203 4 : delete fgTriggerDigits;
204 2 : fgTriggerDigits = 0;
205 :
206 4 : delete fgRawUtils;
207 2 : fgRawUtils = 0;
208 4 : delete fgClusterizer;
209 2 : fgClusterizer = 0;
210 :
211 4 : delete fgTriggerProcessor;
212 2 : fgTriggerProcessor = 0;
213 :
214 10 : if(fMatches) { fMatches->Delete(); delete fMatches; fMatches = 0;}
215 :
216 4 : AliCodeTimer::Instance()->Print();
217 6 : }
218 :
219 : ///
220 : /// Init the fgClusterizer with geometry and calibration pointers, avoid doing it twice.
221 : ///
222 : //____________________________________________________________________________
223 : void AliEMCALReconstructor::InitClusterizer() const
224 : {
225 : Int_t clusterizerType = -1;
226 : Int_t eventType = -1;
227 16 : if(GetRecParam())
228 : {
229 8 : clusterizerType = GetRecParam()->GetClusterizerFlag();
230 8 : eventType = GetRecParam()->GetEventSpecie();
231 8 : }
232 : else
233 : {
234 0 : AliCDBEntry *entry = (AliCDBEntry*)
235 0 : AliCDBManager::Instance()->Get("EMCAL/Calib/RecoParam");
236 :
237 : // Get The reco param for the default event specie
238 0 : if (entry)
239 : {
240 0 : AliEMCALRecParam *recParam = (AliEMCALRecParam*)((TObjArray *) entry->GetObject())->At(0);
241 0 : if(recParam) clusterizerType = recParam->GetClusterizerFlag();
242 0 : }
243 : }
244 :
245 : // Check if clusterizer previously set corresponds to what is needed for this event type
246 8 : if(fgClusterizer)
247 : {
248 6 : if(eventType!=AliRecoParam::kCalib)
249 : {
250 : //printf("ReCreate clusterizer? Clusterizer set <%d>, Clusterizer in use <%s>\n",
251 : // clusterizerType, fgClusterizer->Version());
252 :
253 18 : if (clusterizerType == AliEMCALRecParam::kClusterizerv1 && !strcmp(fgClusterizer->Version(),"clu-v1")) return;
254 :
255 0 : else if(clusterizerType == AliEMCALRecParam::kClusterizerNxN && !strcmp(fgClusterizer->Version(),"clu-NxN")) return;
256 :
257 0 : else if(clusterizerType == AliEMCALRecParam::kClusterizerv2 && !strcmp(fgClusterizer->Version(),"clu-v2")) return;
258 :
259 : //Need to create new clusterizer, the one set previously is not the correct one
260 0 : delete fgClusterizer;
261 : }
262 0 : else return;
263 : }
264 :
265 2 : if (clusterizerType == AliEMCALRecParam::kClusterizerv1)
266 : {
267 4 : fgClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData,fCalibTime,fPedestalData);
268 2 : }
269 0 : else if (clusterizerType == AliEMCALRecParam::kClusterizerNxN)
270 : {
271 0 : fgClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData,fCalibTime,fPedestalData);
272 0 : }
273 0 : else if (clusterizerType == AliEMCALRecParam::kClusterizerv2)
274 : {
275 0 : fgClusterizer = new AliEMCALClusterizerv2 (fGeom, fCalibData,fCalibTime,fPedestalData);
276 0 : }
277 : else
278 : {
279 0 : AliFatal(Form("Unknown clusterizer %d ", clusterizerType));
280 : }
281 10 : }
282 :
283 : ///
284 : /// Method called by AliReconstruction;
285 : /// Only the clusterization is performed, the rest of the reconstruction is done
286 : /// in FillESD because the track matching needs access to the AliESD object to
287 : /// retrieve the tracks reconstructed by the global tracking.
288 : /// Works on the current event.
289 : ///
290 : //____________________________________________________________________________
291 : void AliEMCALReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
292 : {
293 16 : AliCodeTimerAuto("",0)
294 :
295 : // Get input digits and put them in fgDigitsArr, clear the list before
296 8 : ReadDigitsArrayFromTree(digitsTree);
297 :
298 8 : InitClusterizer();
299 :
300 8 : fgClusterizer->InitParameters();
301 8 : fgClusterizer->SetOutput(clustersTree);
302 :
303 : // Skip clusterization of LED events
304 24 : if (GetRecParam()->GetEventSpecie()!=AliRecoParam::kCalib)
305 : {
306 24 : if(fgDigitsArr && fgDigitsArr->GetEntries())
307 : {
308 8 : fgClusterizer->SetInput(digitsTree);
309 :
310 : //fgClusterizer->Digits2Clusters("deb all") ; //For debugging
311 8 : fgClusterizer->Digits2Clusters("");
312 :
313 8 : fgClusterizer->Clear();
314 :
315 : }//digits array exists and has somethind
316 : }//not a LED event
317 :
318 8 : clustersTree->Fill();
319 :
320 : // Deleting the recpoints at the end of the reconstruction call
321 8 : fgClusterizer->DeleteRecPoints();
322 8 : }
323 :
324 : ///
325 : /// Conversion from raw data to EMCAL digits.
326 : /// Works on a single-event basis.
327 : ///
328 : //____________________________________________________________________________
329 : void AliEMCALReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
330 : {
331 8 : rawReader->Reset() ;
332 :
333 16 : for (int i=0;i<fgTriggerData->GetEntriesFast();i++)
334 : {
335 4 : ((AliEMCALTriggerData*)fgTriggerData->At(i))->SetMode(1);
336 : }
337 :
338 8 : if(fgDigitsArr) fgDigitsArr->Clear("C");
339 :
340 4 : const int kNTRU = fGeom->GetNTotalTRU();
341 8 : TClonesArray *digitsTrg = new TClonesArray("AliEMCALTriggerRawDigit", kNTRU * 96);
342 :
343 : Int_t bufsize = 32000;
344 4 : digitsTree->Branch("EMCAL", &fgDigitsArr, bufsize);
345 4 : digitsTree->Branch("EMTRG", &digitsTrg, bufsize);
346 :
347 : // Skip calibration events do the rest
348 : Bool_t doFit = kTRUE;
349 8 : if ( !(GetRecParam()->FitLEDEvents()) && GetRecParam()->GetEventSpecie()==AliRecoParam::kCalib) doFit = kFALSE;
350 :
351 4 : if (doFit)
352 : {
353 : // must be done here because, in constructor, option is not yet known
354 4 : fgRawUtils->SetOption(GetOption());
355 :
356 : // fgRawUtils->SetRawFormatHighLowGainFactor(GetRecParam()->GetHighLowGainFactor());
357 :
358 : // fgRawUtils->SetRawFormatOrder(GetRecParam()->GetOrderParameter());
359 : // fgRawUtils->SetRawFormatTau(GetRecParam()->GetTau());
360 4 : fgRawUtils->SetNoiseThreshold(GetRecParam()->GetNoiseThreshold());
361 4 : fgRawUtils->SetNPedSamples(GetRecParam()->GetNPedSamples());
362 4 : fgRawUtils->SetRemoveBadChannels(GetRecParam()->GetRemoveBadChannels());
363 8 : if (!fgRawUtils->GetFittingAlgorithm()) fgRawUtils->SetFittingAlgorithm(GetRecParam()->GetFittingAlgorithm());
364 4 : fgRawUtils->SetFALTROUsage(GetRecParam()->UseFALTRO());
365 4 : fgRawUtils->SetL1PhaseUsage(GetRecParam()->UseL1Phase());
366 : // fgRawUtils->SetFALTROUsage(0);
367 :
368 : //fgRawUtils->SetTimeMin(GetRecParam()->GetTimeMin());
369 : //fgRawUtils->SetTimeMax(GetRecParam()->GetTimeMax());
370 :
371 : // fgRawUtils->SetTimeMin(-99999 );
372 : // fgRawUtils->SetTimeMax( 99999 );
373 :
374 4 : fgRawUtils->Raw2Digits(rawReader,fgDigitsArr,fPedestalData,digitsTrg,fgTriggerData);
375 :
376 4 : }//skip calibration event
377 : else
378 : {
379 0 : AliDebug(1," Calibration Event, skip!");
380 : }
381 :
382 4 : digitsTree->Fill();
383 4 : digitsTrg->Delete();
384 8 : delete digitsTrg;
385 :
386 4 : }
387 :
388 : ///
389 : /// Create AliESDCaloCells, AliESDCaloClusters and AliESDCaloTrigger entries
390 : /// from AliEMCALRecPoints, AliEMCALDigits.
391 : /// Called by AliReconstruct after Reconstruct() and global tracking and
392 : /// vertexing and V0.
393 : /// Works on the current event
394 : ///
395 : //____________________________________________________________________________
396 : void AliEMCALReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
397 : AliESDEvent* esd) const
398 : {
399 : //########################################
400 : // Trigger
401 : //########################################
402 :
403 : static int saveOnce[2] = {0,0};
404 :
405 16 : Int_t v0M[2] = {0, 0};
406 :
407 8 : AliESDVZERO* esdV0 = esd->GetVZEROData();
408 :
409 8 : if (esdV0)
410 : {
411 8 : v0M[0] = esdV0->GetTriggerChargeC();
412 8 : v0M[1] = esdV0->GetTriggerChargeA();
413 8 : }
414 : else
415 : {
416 0 : AliWarning("No V0 ESD! Run trigger processor w/ null V0 charges");
417 : }
418 :
419 19 : if (fgTriggerDigits && fgTriggerDigits->GetEntriesFast()) fgTriggerDigits->Delete();
420 :
421 8 : TBranch *branchtrg = digitsTree->GetBranch("EMTRG");
422 :
423 8 : if (!branchtrg)
424 : {
425 0 : AliError("Can't get the branch with the EMCAL trigger digits!");
426 0 : return;
427 : }
428 :
429 8 : branchtrg->SetAddress(&fgTriggerDigits);
430 8 : branchtrg->GetEntry(0);
431 :
432 : // Note: fgTriggerProcessor reset done at the end of this method
433 :
434 : // TO DO: Run2 emulation
435 8 : if(esd->GetRunNumber() < 200000) // Enable trigger emulation for Run1
436 8 : fgTriggerProcessor->Digits2Trigger(fgTriggerDigits, v0M, (AliEMCALTriggerData*)fgTriggerData->At(0));
437 :
438 : // Fill ESD
439 16 : AliESDCaloTrigger* trgESD = esd->GetCaloTrigger("EMCAL");
440 :
441 8 : if (trgESD)
442 : {
443 8 : trgESD->Allocate(fgTriggerDigits->GetEntriesFast());
444 : int trgBitWord=0;
445 300 : for (Int_t i = 0; i < fgTriggerDigits->GetEntriesFast(); i++) {
446 142 : AliEMCALTriggerRawDigit* rdig = (AliEMCALTriggerRawDigit*)fgTriggerDigits->At(i);
447 426 : if (AliDebugLevel() > 999) rdig->Print("");
448 :
449 142 : Int_t px, py, trgBitWord=0;
450 142 : if (fGeom->GetPositionInEMCALFromAbsFastORIndex(rdig->GetId(), px, py)) {
451 142 : Int_t a = -1, t = -1, times[10];
452 :
453 142 : rdig->GetMaximum(a, t);
454 142 : rdig->GetL0Times(times);
455 :
456 142 : trgESD->Add(px, py, a, t, times, rdig->GetNL0Times(), rdig->GetL1TimeSum(), rdig->GetL1SubRegion(), rdig->GetTriggerBits());
457 142 : trgBitWord |= rdig->GetTriggerBits();
458 168 : if (rdig->GetNL0Times()) trgBitWord |= 1 << 5; // Overwrite L0 id from STU payload
459 142 : }
460 142 : }
461 :
462 8 : trgESD->SetTriggerBitWord(trgBitWord);
463 :
464 40 : for (int i=0;i<fgTriggerData->GetEntriesFast();i++)
465 : {
466 72 : for (int j=0;j<2;j++) {
467 24 : trgESD->SetL1Threshold(i,2 * j , ((AliEMCALTriggerData*)fgTriggerData->At(i))->GetL1JetThreshold( j));
468 24 : trgESD->SetL1Threshold(i,2 * j + 1, ((AliEMCALTriggerData*)fgTriggerData->At(i))->GetL1GammaThreshold(j));
469 : }
470 12 : trgESD->SetL1FrameMask(i,((AliEMCALTriggerData*)fgTriggerData->At(i))->GetL1FrameMask());
471 :
472 12 : Int_t v0[2];
473 12 : ((AliEMCALTriggerData*)fgTriggerData->At(i))->GetL1V0(v0);
474 12 : trgESD->SetL1V0(i,v0);
475 :
476 12 : trgESD->SetMedian(i,((AliEMCALTriggerData*)fgTriggerData->At(i))->GetMedian());
477 24 : if (!saveOnce[i] && ((AliEMCALTriggerData*)fgTriggerData->At(i))->GetL1DataDecoded())
478 : {
479 0 : int type[19] = {0};
480 0 : ((AliEMCALTriggerData*)fgTriggerData->At(i))->GetL1TriggerType(type);
481 :
482 0 : esd->SetCaloTriggerType(i,type);
483 :
484 0 : saveOnce[i] = 1;
485 0 : }
486 12 : }
487 8 : }
488 :
489 : // Resetting
490 40 : for (int i=0;i<fgTriggerData->GetEntriesFast();i++)
491 : {
492 12 : ((AliEMCALTriggerData*)fgTriggerData->At(i))->Reset();
493 : }
494 :
495 : //########################################
496 : //##############Fill CaloCells############
497 : //########################################
498 :
499 : // Get input digits and put them in fgDigitsArr, clear the list before
500 8 : ReadDigitsArrayFromTree(digitsTree);
501 :
502 8 : Int_t nDigits = fgDigitsArr->GetEntries(), idignew = 0 ;
503 24 : AliDebug(1,Form("%d digits",nDigits));
504 :
505 8 : AliESDCaloCells &emcCells = *(esd->GetEMCALCells());
506 8 : emcCells.CreateContainer(nDigits);
507 8 : emcCells.SetType(AliVCaloCells::kEMCALCell);
508 :
509 8 : Float_t energy = 0;
510 8 : Float_t time = 0;
511 8 : Int_t mapDigitAndCellIndex[20000] = { 0 } ; // needed to pack cell MC labels in cluster and not loop all digits
512 :
513 246 : for (Int_t idig = 0 ; idig < nDigits ; idig++)
514 : {
515 115 : const AliEMCALDigit * dig = (const AliEMCALDigit*)fgDigitsArr->At(idig);
516 115 : time = dig->GetTime(); // Time already calibrated in clusterizer
517 115 : energy = dig->GetAmplitude(); // energy calibrated in clusterizer
518 :
519 115 : if (energy <= 0 ) continue ;
520 :
521 115 : fgClusterizer->Calibrate(energy,time,dig->GetId()); // Digits already calibrated in clusterizers
522 :
523 117 : if (energy <= 0 ) continue ; // Digits tagged as bad (dead, hot, not alive) are set to 0 in calibrate, remove them
524 :
525 : // Only for MC
526 : // Get the label of the primary particle that generated the cell
527 : // Assign the particle that deposited more energy
528 113 : Int_t nparents = dig->GetNiparent();
529 : Int_t digLabel =-1 ;
530 : Float_t edep =-1.;
531 113 : if ( nparents > 0 )
532 : {
533 202 : for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
534 : {
535 51 : if(edep >= dig->GetDEParent(jndex+1)) continue ;
536 :
537 50 : digLabel = dig->GetIparent (jndex+1);
538 50 : edep = dig->GetDEParent(jndex+1);
539 50 : } // all primaries in digit
540 50 : } // select primary label
541 :
542 : Bool_t highGain = kFALSE;
543 226 : if( dig->GetType() == AliEMCALDigit::kHG ) highGain = kTRUE;
544 :
545 113 : emcCells.SetCell(idignew,dig->GetId(),energy, time,digLabel,0.,highGain);
546 :
547 113 : mapDigitAndCellIndex[idignew] = idig; // needed to pack cell MC labels in cluster
548 :
549 113 : idignew++;
550 113 : }
551 :
552 8 : emcCells.SetNumberOfCells(idignew);
553 8 : emcCells.Sort();
554 :
555 : //------------------------------------------------------------
556 : //-----------------CLUSTERS-----------------------------------
557 : //------------------------------------------------------------
558 8 : clustersTree->SetBranchStatus("*",0); //disable all branches
559 8 : clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
560 16 : if(fgClustersArr) fgClustersArr->Clear();
561 8 : TBranch *branch = clustersTree->GetBranch("EMCALECARP");
562 8 : branch->SetAddress(&fgClustersArr);
563 8 : branch->GetEntry(0);
564 : //clustersTree->GetEvent(0);
565 :
566 8 : Int_t nClusters = fgClustersArr->GetEntries(), nClustersNew=0;
567 24 : AliDebug(1,Form("%d clusters",nClusters));
568 :
569 :
570 : //########################################
571 : //##############Fill CaloClusters#########
572 : //########################################
573 82 : for (Int_t iClust = 0 ; iClust < nClusters ; iClust++)
574 : {
575 33 : const AliEMCALRecPoint * clust = (const AliEMCALRecPoint*)fgClustersArr->At(iClust);
576 :
577 33 : if(!clust) continue;
578 :
579 : // clust->Print(); //For debugging
580 :
581 : // Get information from EMCAL reconstruction points
582 33 : Float_t xyz[3];
583 33 : TVector3 gpos;
584 33 : clust->GetGlobalPosition(gpos);
585 264 : for (Int_t ixyz=0; ixyz<3; ixyz++)
586 198 : xyz[ixyz] = gpos[ixyz];
587 33 : Float_t elipAxis[2];
588 33 : clust->GetElipsAxis(elipAxis);
589 :
590 : //Create digits lists
591 33 : Int_t cellMult = clust->GetMultiplicity();
592 :
593 33 : Float_t *amplFloat = clust->GetEnergiesList();
594 33 : Int_t *digitInts = clust->GetAbsId();
595 33 : TArrayS absIdList(cellMult);
596 33 : TArrayD fracList(cellMult);
597 :
598 : Int_t newCellMult = 0;
599 260 : for (Int_t iCell=0; iCell<cellMult; iCell++)
600 : {
601 97 : if (amplFloat[iCell] > 0)
602 : {
603 194 : absIdList[newCellMult] = (UShort_t)(digitInts[iCell]);
604 :
605 : //Calculate Fraction
606 388 : if(emcCells.GetCellAmplitude(digitInts[iCell])>0 && GetRecParam()->GetUnfold())
607 0 : fracList[newCellMult] = amplFloat[iCell]/(emcCells.GetCellAmplitude(digitInts[iCell]));//get cell calibration value
608 : else
609 194 : fracList[newCellMult] = 0;
610 :
611 97 : newCellMult++;
612 97 : }
613 : }
614 :
615 33 : absIdList.Set(newCellMult);
616 33 : fracList.Set(newCellMult);
617 :
618 : // Accept cluster if it has some digit
619 33 : if(newCellMult > 0)
620 : {
621 33 : nClustersNew++;
622 :
623 : // Fill the ESDCaloCluster
624 66 : AliESDCaloCluster * ec = new AliESDCaloCluster() ;
625 33 : ec->SetType(AliVCluster::kEMCALClusterv1);
626 33 : ec->SetPosition(xyz);
627 66 : ec->SetE(clust->GetEnergy());
628 33 : ec->SetTOF(clust->GetTime()) ; //time-of-flight
629 :
630 33 : ec->SetChi2(-1); //not yet implemented
631 :
632 : // Distance to the nearest bad tower
633 33 : ec->SetDistanceToBadChannel(clust->GetDistanceToBadTower());
634 :
635 : //
636 : // Shower shape related
637 : //
638 33 : ec->SetNCells(newCellMult);
639 66 : ec->SetDispersion(clust->GetDispersion());
640 33 : ec->SetM02(elipAxis[0]*elipAxis[0]) ;
641 33 : ec->SetM20(elipAxis[1]*elipAxis[1]) ;
642 33 : ec->SetNExMax(clust->GetNExMax()); //number of local maxima
643 :
644 : //
645 : // Primaries
646 : //
647 33 : Int_t parentMult = 0;
648 33 : Int_t *parentList = clust->GetParents(parentMult); // label index
649 33 : Float_t *parentListDE = clust->GetParentsDE(); // deposited energy
650 :
651 33 : ec->SetLabel(parentList,parentMult);
652 :
653 33 : ec->SetClusterMCEdepFractionFromEdepArray(parentListDE);
654 :
655 : //
656 : // Change type of cells list from short to ushort
657 : // Pack the at maximum 4 different primary contributions on cells into a single int.
658 : // (Profit from the same loop on cells in cluster)
659 : //
660 66 : UShort_t *newAbsIdList = new UShort_t[newCellMult];
661 66 : Double_t *newFracList = new Double_t[newCellMult];
662 66 : UInt_t *cellEDepFrac = new UInt_t [newCellMult];
663 :
664 260 : for(Int_t i = 0; i < newCellMult ; i++)
665 : {
666 : // int to short
667 194 : newAbsIdList[i] = absIdList[i];
668 194 : newFracList [i] = fracList[i];
669 97 : cellEDepFrac[i] = 0;
670 :
671 : // *Pack cells parent energy deposit*
672 97 : if(parentMult > 0)
673 : {
674 : // Get the digit that originated this cell cluster
675 98 : Int_t cellPos = emcCells.GetCellPosition (newAbsIdList[i]);
676 49 : Int_t idigit = mapDigitAndCellIndex[cellPos];
677 98 : const AliEMCALDigit * dig = (const AliEMCALDigit*)fgDigitsArr->At(idigit);
678 :
679 : // Find the 4 MC labels that contributed to the cluster and their
680 : // deposited energy in the current digit
681 49 : Int_t nparents = dig->GetNiparent();
682 49 : if ( nparents > 0 )
683 : {
684 : Int_t digLabel =-1 ;
685 : Float_t edep = 0 ;
686 : Float_t edepTot = 0 ;
687 44 : Float_t mcEDepFrac[4] = {0,0,0,0};
688 :
689 178 : for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
690 : { // all primaries in digit
691 45 : digLabel = dig->GetIparent (jndex+1);
692 45 : edep = dig->GetDEParent(jndex+1);
693 45 : edepTot += edep;
694 :
695 89 : if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
696 2 : else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
697 0 : else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
698 0 : else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
699 : } // all primaries in digit
700 :
701 : // Divide energy deposit by total deposited energy
702 : // Do not take the stored digit energy since it is smeared and with noise
703 : // One could go back to the sdigit, but it is simpler to take the added
704 : // deposited energy of all primaries.
705 : // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
706 44 : if(edepTot > 0.01)
707 : {
708 440 : for(Int_t idep = 0; idep < 4; idep++)
709 176 : mcEDepFrac[idep] /= edepTot;
710 :
711 88 : cellEDepFrac[i] = ec->PackMCEdepFraction(mcEDepFrac);
712 44 : }
713 44 : } // at least one parent in digit
714 49 : } // at least one parent in cluster, do the cell primary packing
715 : } // cell cluster loop
716 :
717 33 : ec->SetCellsAbsId(newAbsIdList);
718 33 : ec->SetCellsAmplitudeFraction(newFracList);
719 33 : ec->SetCellsMCEdepFractionMap(cellEDepFrac);
720 :
721 : //
722 : // Track matching
723 : //
724 33 : fMatches->Clear();
725 33 : Int_t nTracks = esd->GetNumberOfTracks();
726 1276 : for (Int_t itrack = 0; itrack < nTracks; itrack++)
727 : {
728 605 : AliESDtrack * track = esd->GetTrack(itrack) ; // retrieve track
729 1210 : if(track->GetEMCALcluster()==iClust)
730 : {
731 18 : Float_t dEta=-999, dPhi=-999;
732 36 : Bool_t isMatch = CalculateResidual(track, ec, dEta, dPhi);
733 18 : if(!isMatch)
734 : {
735 : // AliDebug(10, "Not good");
736 0 : continue;
737 : }
738 :
739 36 : AliEMCALMatch *match = new AliEMCALMatch();
740 18 : match->SetIndexT(itrack);
741 18 : match->SetDistance(TMath::Sqrt(dEta*dEta+dPhi*dPhi));
742 18 : match->SetdEta(dEta);
743 18 : match->SetdPhi(dPhi);
744 18 : fMatches->Add(match);
745 36 : }
746 605 : }
747 :
748 33 : fMatches->Sort(kSortAscending); //Sort matched tracks from closest to furthest
749 33 : Int_t nMatch = fMatches->GetEntries();
750 33 : TArrayI arrayTrackMatched(nMatch);
751 102 : for(Int_t imatch=0; imatch<nMatch; imatch++)
752 : {
753 36 : AliEMCALMatch *match = (AliEMCALMatch*)fMatches->At(imatch);
754 36 : arrayTrackMatched[imatch] = match->GetIndexT();
755 18 : if(imatch==0)
756 : {
757 14 : ec->SetTrackDistance(match->GetdPhi(), match->GetdEta());
758 : }
759 : }
760 :
761 33 : ec->AddTracksMatched(arrayTrackMatched);
762 :
763 : //add the cluster to the esd object
764 33 : esd->AddCaloCluster(ec);
765 :
766 66 : delete ec;
767 66 : delete [] newAbsIdList ;
768 66 : delete [] newFracList ;
769 66 : delete [] cellEDepFrac ;
770 33 : }
771 33 : } // cycle on clusters
772 :
773 : //
774 : // Reset the index of matched cluster for tracks
775 : // to the one in CaloCluster array
776 : //
777 8 : Int_t ncls = esd->GetNumberOfCaloClusters();
778 102 : for(Int_t icl=0; icl<ncls; icl++)
779 : {
780 43 : AliESDCaloCluster *cluster = esd->GetCaloCluster(icl);
781 96 : if(!cluster || !cluster->IsEMCAL()) continue;
782 33 : TArrayI *trackIndex = cluster->GetTracksMatched();
783 102 : for(Int_t itr=0; itr<trackIndex->GetSize(); itr++)
784 : {
785 18 : AliESDtrack *track = esd->GetTrack(trackIndex->At(itr));
786 18 : track->SetEMCALcluster(cluster->GetID());
787 : }
788 33 : }
789 :
790 : // Fill ESDCaloCluster with PID weights
791 8 : AliEMCALPID *pid = new AliEMCALPID;
792 : //pid->SetPrintInfo(kTRUE);
793 8 : pid->SetReconstructor(kTRUE);
794 8 : pid->RunPID(esd);
795 16 : delete pid;
796 :
797 : // Store EMCAL misalignment matrixes
798 8 : FillMisalMatrixes(esd) ;
799 :
800 16 : }
801 :
802 : ///
803 : /// Store EMCAL matrixes in ESD Header.
804 : ///
805 : //==================================================================================
806 : void AliEMCALReconstructor::FillMisalMatrixes(AliESDEvent* esd)const
807 : {
808 : // Check, if matrixes was already stored
809 86 : for(Int_t sm = 0 ; sm < fGeom->GetNumberOfSuperModules(); sm++)
810 : {
811 36 : if(esd->GetEMCALMatrix(sm)!=0)
812 6 : return ;
813 : }
814 :
815 : // Create and store matrixes
816 2 : if(!gGeoManager)
817 : {
818 0 : AliError("Can not store misal. matrixes: no gGeoManager! \n") ;
819 0 : return ;
820 : }
821 :
822 : // Note, that owner of copied matrixes will be header
823 : const Int_t bufsize = 255;
824 2 : char path[bufsize] ;
825 : TGeoHMatrix * m = 0x0;
826 : Int_t tmpType = -1;
827 : Int_t SMOrder = 0;
828 2 : TString SMName;
829 96 : for(Int_t sm = 0; sm < fGeom->GetNumberOfSuperModules(); sm++)
830 : {
831 80 : if(fGeom->GetSMType(sm) == AliEMCALGeometry::kEMCAL_Standard ) SMName = "SMOD";
832 20 : else if(fGeom->GetSMType(sm) == AliEMCALGeometry::kEMCAL_Half ) SMName = "SM10";
833 22 : else if(fGeom->GetSMType(sm) == AliEMCALGeometry::kEMCAL_3rd ) SMName = "SM3rd";
834 22 : else if( fGeom->GetSMType(sm) == AliEMCALGeometry::kDCAL_Standard ) SMName = "DCSM";
835 6 : else if( fGeom->GetSMType(sm) == AliEMCALGeometry::kDCAL_Ext ) SMName = "DCEXT";
836 0 : else AliError("Unkown SM Type!!");
837 :
838 60 : if(fGeom->GetSMType(sm) == tmpType)
839 : {
840 25 : SMOrder++;
841 25 : }
842 : else
843 : {
844 5 : tmpType = fGeom->GetSMType(sm);
845 : SMOrder = 1;
846 : }
847 60 : snprintf(path,bufsize,"/ALIC_1/XEN1_1/%s_%d", SMName.Data(), SMOrder) ;
848 :
849 60 : if (gGeoManager->CheckPath(path))
850 : {
851 30 : gGeoManager->cd(path);
852 30 : m = gGeoManager->GetCurrentMatrix() ;
853 : // printf("================================================= \n");
854 : // printf("AliEMCALReconstructor::FixMisalMatrixes(), sm %d, \n",sm);
855 : // m->Print("");
856 90 : esd->SetEMCALMatrix(new TGeoHMatrix(*m),sm) ;
857 : // printf("================================================= \n");
858 : }
859 : else
860 : {
861 0 : esd->SetEMCALMatrix(NULL,sm) ;
862 : }
863 : }
864 10 : }
865 :
866 : ///
867 : /// Read the digits from the input tree
868 : /// See AliEMCALClusterizer::SetInput(TTree *digitsTree);
869 : ///
870 : //__________________________________________________________________________
871 : void AliEMCALReconstructor::ReadDigitsArrayFromTree(TTree *digitsTree) const
872 : {
873 : // Clear previous digits in the list
874 32 : if(fgDigitsArr)
875 : {
876 16 : fgDigitsArr->Clear("C");
877 16 : }
878 : else
879 : {
880 : // It should not happen, but just in case ...
881 0 : fgDigitsArr = new TClonesArray("AliEMCALDigit",100);
882 : }
883 :
884 : // Read the digits from the input tree
885 16 : TBranch *branch = digitsTree->GetBranch("EMCAL");
886 16 : if (!branch)
887 : {
888 0 : AliError("can't get the branch with the EMCAL digits !");
889 0 : return;
890 : }
891 :
892 16 : branch->SetAddress(&fgDigitsArr);
893 16 : branch->GetEntry(0);
894 32 : }
895 :
896 : ///
897 : /// Calculate the residual between track and cluster.
898 : ///
899 : /// If the esdFriend is available, use the TPCOuter point as the starting
900 : /// point of extrapolation, otherwise use the TPCInner point.
901 : ///
902 : //==================================================================================
903 : Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCluster *cluster, Float_t &dEta, Float_t &dPhi)const
904 : {
905 36 : dEta = -999, dPhi = -999;
906 : Bool_t ITSTrackSA = 0;
907 :
908 : AliExternalTrackParam *trkParam = 0;
909 :
910 18 : const AliESDfriendTrack* friendTrack = track->GetFriendTrack();
911 36 : if(friendTrack && friendTrack->GetTPCOut())
912 18 : trkParam = const_cast<AliExternalTrackParam*>(friendTrack->GetTPCOut());
913 0 : else if(track->GetInnerParam())
914 0 : trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
915 : else
916 : {
917 0 : trkParam = new AliExternalTrackParam(*track); //If there is ITSSa track
918 : ITSTrackSA = 1;
919 : }
920 18 : if(!trkParam) return kFALSE;
921 :
922 18 : AliExternalTrackParam trkParamTmp (*trkParam);
923 72 : if(!AliEMCALRecoUtils::ExtrapolateTrackToCluster(&trkParamTmp, cluster, track->GetMass(kTRUE), GetRecParam()->GetExtrapolateStep(), dEta, dPhi))
924 : {
925 0 : if(ITSTrackSA) delete trkParam;
926 0 : return kFALSE;
927 : }
928 :
929 36 : if(ITSTrackSA) delete trkParam;
930 18 : return kTRUE;
931 36 : }
932 :
933 : ///
934 : /// Default constructor.
935 : ///
936 : //==================================================================================
937 : AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch()
938 18 : : TObject(),
939 18 : fIndexT(-1),
940 18 : fDistance(-999.),
941 18 : fdEta(-999.),
942 18 : fdPhi(-999.)
943 90 : {
944 36 : }
945 :
946 : ///
947 : /// Copy constructor.
948 : ///
949 : //==================================================================================
950 : AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
951 0 : : TObject(),
952 0 : fIndexT(copy.fIndexT),
953 0 : fDistance(copy.fDistance),
954 0 : fdEta(copy.fdEta),
955 0 : fdPhi(copy.fdPhi)
956 0 : {
957 0 : }
958 :
959 : ///
960 : /// Assignment operator; use copy ctor.
961 : ///
962 : //_____________________________________________________________________
963 : AliEMCALReconstructor::AliEMCALMatch& AliEMCALReconstructor::AliEMCALMatch::AliEMCALMatch::operator = (const AliEMCALMatch &source)
964 : {
965 0 : if (&source == this) return *this;
966 :
967 0 : new (this) AliEMCALMatch(source);
968 0 : return *this;
969 0 : }
970 :
971 : ///
972 : /// Compare wrt the residual.
973 : ///
974 : //==================================================================================
975 : Int_t AliEMCALReconstructor::AliEMCALMatch::Compare(const TObject *obj) const
976 : {
977 8 : AliEMCALReconstructor::AliEMCALMatch *that = (AliEMCALReconstructor::AliEMCALMatch*)obj;
978 :
979 4 : Double_t thisDist = fDistance;//fDistance;
980 4 : Double_t thatDist = that->fDistance;//that->GetDistance();
981 :
982 8 : if (thisDist > thatDist) return 1;
983 0 : else if (thisDist < thatDist) return -1;
984 0 : return 0;
985 4 : }
|
