]>
Commit | Line | Data |
---|---|---|
6fdd30c4 | 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 | /* $Id: $ */ | |
16 | ||
17 | //_________________________________________________________________________ | |
18 | // This analysis provides a new list of clusters to be used in other analysis | |
19 | // | |
20 | // Author: Gustavo Conesa Balbastre, | |
21 | // Adapted from analysis class from Deepa Thomas | |
22 | // | |
23 | // | |
24 | //_________________________________________________________________________ | |
25 | ||
26 | // --- Root --- | |
27 | #include "TString.h" | |
28 | #include "TRefArray.h" | |
29 | #include "TClonesArray.h" | |
30 | #include "TTree.h" | |
31 | #include "TGeoManager.h" | |
5994e71f | 32 | #include "TROOT.h" |
33 | #include "TInterpreter.h" | |
6fdd30c4 | 34 | |
35 | // --- AliRoot Analysis Steering | |
36 | #include "AliAnalysisTask.h" | |
37 | #include "AliAnalysisManager.h" | |
38 | #include "AliESDEvent.h" | |
39 | #include "AliGeomManager.h" | |
40 | #include "AliVCaloCells.h" | |
41 | #include "AliAODCaloCluster.h" | |
42 | #include "AliCDBManager.h" | |
43 | #include "AliCDBStorage.h" | |
44 | #include "AliCDBEntry.h" | |
45 | #include "AliLog.h" | |
46 | #include "AliVEventHandler.h" | |
6060ed91 | 47 | #include "AliAODInputHandler.h" |
6fdd30c4 | 48 | |
49 | // --- EMCAL | |
50 | #include "AliEMCALRecParam.h" | |
51 | #include "AliEMCALAfterBurnerUF.h" | |
52 | #include "AliEMCALGeometry.h" | |
53 | #include "AliEMCALClusterizerNxN.h" | |
54 | #include "AliEMCALClusterizerv1.h" | |
55 | #include "AliEMCALRecPoint.h" | |
56 | #include "AliEMCALDigit.h" | |
57 | #include "AliCaloCalibPedestal.h" | |
58 | #include "AliEMCALCalibData.h" | |
385b7abf | 59 | #include "AliEMCALRecoUtils.h" |
6fdd30c4 | 60 | |
61 | #include "AliAnalysisTaskEMCALClusterize.h" | |
62 | ||
63 | ClassImp(AliAnalysisTaskEMCALClusterize) | |
64 | ||
65 | //________________________________________________________________________ | |
66 | AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) | |
67 | : AliAnalysisTaskSE(name) | |
68 | , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE) | |
69 | , fCalibData(0), fPedestalData(0), fOCDBpath("raw://") | |
70 | , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0) | |
71 | , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE) | |
72 | , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kTRUE) | |
6060ed91 | 73 | , fRun(-1), fRecoUtils(0), fConfigName(""), fCellLabels() |
6fdd30c4 | 74 | |
75 | { | |
76 | //ctor | |
6060ed91 | 77 | for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0; |
78 | for(Int_t j = 0; j < 12672; j++) fCellLabels[j] = -1; | |
6fdd30c4 | 79 | fDigitsArr = new TClonesArray("AliEMCALDigit",200); |
80 | fClusterArr = new TObjArray(100); | |
81 | fCaloClusterArr = new TObjArray(100); | |
82 | fRecParam = new AliEMCALRecParam; | |
385b7abf | 83 | fBranchNames = "ESD:AliESDHeader.,EMCALCells."; |
84 | fRecoUtils = new AliEMCALRecoUtils(); | |
6fdd30c4 | 85 | } |
86 | ||
87 | //________________________________________________________________________ | |
88 | AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() | |
89 | : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize") | |
90 | , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE) | |
91 | , fCalibData(0), fPedestalData(0), fOCDBpath("raw://") | |
92 | , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0) | |
542e1379 | 93 | , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE) |
6fdd30c4 | 94 | , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kFALSE) |
6060ed91 | 95 | , fRun(-1), fRecoUtils(0), fConfigName(""), fCellLabels() |
6fdd30c4 | 96 | { |
97 | // Constructor | |
6060ed91 | 98 | for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0; |
99 | for(Int_t j = 0; j < 24*48*11; j++) fCellLabels[j] = -1; | |
6fdd30c4 | 100 | fDigitsArr = new TClonesArray("AliEMCALDigit",200); |
101 | fClusterArr = new TObjArray(100); | |
102 | fCaloClusterArr = new TObjArray(100); | |
103 | fRecParam = new AliEMCALRecParam; | |
385b7abf | 104 | fBranchNames = "ESD:AliESDHeader.,EMCALCells."; |
105 | fRecoUtils = new AliEMCALRecoUtils(); | |
6fdd30c4 | 106 | } |
107 | ||
5994e71f | 108 | |
6fdd30c4 | 109 | //________________________________________________________________________ |
110 | AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize() | |
111 | { | |
112 | //dtor | |
113 | ||
114 | if (fDigitsArr){ | |
115 | fDigitsArr->Clear("C"); | |
542e1379 | 116 | delete fDigitsArr; |
6fdd30c4 | 117 | } |
118 | ||
119 | if (fClusterArr){ | |
120 | fClusterArr->Delete(); | |
542e1379 | 121 | delete fClusterArr; |
6fdd30c4 | 122 | } |
123 | ||
124 | if (fCaloClusterArr){ | |
125 | fCaloClusterArr->Delete(); | |
542e1379 | 126 | delete fCaloClusterArr; |
6fdd30c4 | 127 | } |
128 | ||
385b7abf | 129 | if(fClusterizer) delete fClusterizer; |
130 | if(fUnfolder) delete fUnfolder; | |
131 | if(fRecoUtils) delete fRecoUtils; | |
132 | ||
6fdd30c4 | 133 | } |
134 | ||
5994e71f | 135 | //------------------------------------------------------------------- |
136 | void AliAnalysisTaskEMCALClusterize::Init() | |
137 | { | |
138 | //Init analysis with configuration macro if available | |
139 | ||
140 | if(gROOT->LoadMacro(fConfigName) >=0){ | |
141 | printf("Configure analysis with %s\n",fConfigName.Data()); | |
142 | AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()"); | |
143 | fGeomName = clus->fGeomName; | |
144 | fLoadGeomMatrices = clus->fLoadGeomMatrices; | |
145 | fOCDBpath = clus->fOCDBpath; | |
146 | fRecParam = clus->fRecParam; | |
147 | fJustUnfold = clus->fJustUnfold; | |
148 | fFillAODFile = clus->fFillAODFile; | |
149 | fRecoUtils = clus->fRecoUtils; | |
150 | fConfigName = clus->fConfigName; | |
151 | fOutputAODBranchName = clus->fOutputAODBranchName; | |
152 | for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ; | |
153 | ||
154 | } | |
155 | ||
156 | } | |
157 | ||
6fdd30c4 | 158 | //------------------------------------------------------------------- |
159 | void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects() | |
160 | { | |
161 | // Init geometry, create list of output clusters | |
5994e71f | 162 | |
6fdd30c4 | 163 | fGeom = AliEMCALGeometry::GetInstance(fGeomName) ; |
164 | ||
165 | fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0); | |
166 | fOutputAODBranch->SetName(fOutputAODBranchName); | |
167 | AddAODBranch("TClonesArray", &fOutputAODBranch); | |
5994e71f | 168 | |
6fdd30c4 | 169 | } |
170 | ||
171 | //________________________________________________________________________ | |
172 | void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) | |
173 | { | |
174 | // Main loop | |
175 | // Called for each event | |
176 | ||
5994e71f | 177 | //printf("--- Event %d -- \n",(Int_t)Entry()); |
6fdd30c4 | 178 | //Remove the contents of output list set in the previous event |
179 | fOutputAODBranch->Clear("C"); | |
180 | ||
6060ed91 | 181 | //Get the event |
182 | AliVEvent * event = 0; | |
183 | AliESDEvent * esdevent = 0; | |
184 | ||
185 | //Check if input event are embedded events | |
186 | //If so, take output event | |
187 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); | |
188 | if (aodIH && aodIH->GetMergeEvents()) { | |
189 | //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n"); | |
190 | event = AODEvent(); | |
191 | // printf("InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(), | |
192 | // InputEvent()->GetEMCALCells()->GetNumberOfCells()); | |
193 | // printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(), | |
194 | // AODEvent()->GetEMCALCells()->GetNumberOfCells()); | |
195 | } | |
196 | else { | |
197 | event = InputEvent(); | |
198 | esdevent = dynamic_cast<AliESDEvent*>(InputEvent()); | |
199 | } | |
200 | ||
6fdd30c4 | 201 | if (!event) { |
3fa4fb85 | 202 | Error("UserExec","Event not available"); |
6fdd30c4 | 203 | return; |
204 | } | |
205 | ||
206 | //Magic line to write events to AOD file | |
207 | AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile); | |
3fa4fb85 | 208 | LoadBranches(); |
bd9e8ebd | 209 | |
210 | AccessOCDB(); | |
3fa4fb85 | 211 | |
6fdd30c4 | 212 | //------------------------------------------------------------------------------------- |
213 | //Set the geometry matrix, for the first event, skip the rest | |
214 | //------------------------------------------------------------------------------------- | |
215 | if(!fGeomMatrixSet){ | |
216 | if(fLoadGeomMatrices){ | |
217 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ | |
218 | if(fGeomMatrix[mod]){ | |
219 | if(DebugLevel() > 1) | |
220 | fGeomMatrix[mod]->Print(); | |
221 | fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ; | |
222 | } | |
223 | fGeomMatrixSet=kTRUE; | |
224 | }//SM loop | |
225 | }//Load matrices | |
226 | else if(!gGeoManager){ | |
3fa4fb85 | 227 | Info("UserExec","Get geo matrices from data"); |
6fdd30c4 | 228 | //Still not implemented in AOD, just a workaround to be able to work at least with ESDs |
229 | if(!strcmp(event->GetName(),"AliAODEvent")) { | |
230 | if(DebugLevel() > 1) | |
3fa4fb85 | 231 | Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs."); |
6fdd30c4 | 232 | }//AOD |
233 | else { | |
3fa4fb85 | 234 | if(DebugLevel() > 1) |
235 | Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices."); | |
6fdd30c4 | 236 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ; |
237 | if(!esd) { | |
3fa4fb85 | 238 | Error("UserExec","This event does not contain ESDs?"); |
6fdd30c4 | 239 | return; |
240 | } | |
241 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){ | |
242 | if(DebugLevel() > 1) | |
243 | esd->GetEMCALMatrix(mod)->Print(); | |
244 | if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ; | |
245 | } | |
246 | fGeomMatrixSet=kTRUE; | |
247 | }//ESD | |
248 | }//Load matrices from Data | |
5994e71f | 249 | |
250 | //Recover time dependent corrections, put then in recalibration histograms. Do it once | |
251 | fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber()); | |
252 | ||
6fdd30c4 | 253 | }//first event |
254 | ||
255 | //Get the list of cells needed for unfolding and reclustering | |
256 | AliVCaloCells *cells= event->GetEMCALCells(); | |
257 | Int_t nClustersOrg = 0; | |
258 | //------------------------------------------- | |
259 | //---------Unfolding clusters---------------- | |
260 | //------------------------------------------- | |
261 | if (fJustUnfold) { | |
262 | ||
263 | //Fill the array with the EMCAL clusters, copy them | |
264 | for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++) | |
265 | { | |
266 | AliVCluster *clus = event->GetCaloCluster(i); | |
267 | if(clus->IsEMCAL()){ | |
3fa4fb85 | 268 | //printf("Org Cluster %d, E %f\n",i, clus->E()); |
5994e71f | 269 | |
270 | //recalibrate/remove bad channels/etc if requested | |
271 | if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){ | |
272 | //printf("Remove cluster\n"); | |
273 | continue; | |
274 | } | |
275 | ||
276 | if(fRecoUtils->IsRecalibrationOn()){ | |
277 | //printf("Energy before %f ",clus->E()); | |
278 | fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells); | |
279 | //printf("; after %f\n",clus->E()); | |
280 | } | |
281 | //Cast to ESD or AOD, needed to create the cluster array | |
e615d952 | 282 | AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus); |
283 | AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus); | |
284 | if (esdCluster){ | |
285 | fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) ); | |
6fdd30c4 | 286 | }//ESD |
e615d952 | 287 | else if(aodCluster){ |
288 | fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) ); | |
6fdd30c4 | 289 | }//AOD |
3fa4fb85 | 290 | else |
291 | Warning("UserExec()"," - Wrong CaloCluster type?"); | |
6fdd30c4 | 292 | nClustersOrg++; |
293 | } | |
294 | } | |
295 | ||
296 | //Do the unfolding | |
297 | fUnfolder->UnfoldClusters(fCaloClusterArr, cells); | |
298 | ||
299 | //CLEAN-UP | |
300 | fUnfolder->Clear(); | |
301 | ||
5994e71f | 302 | |
6fdd30c4 | 303 | //printf("nClustersOrg %d\n",nClustersOrg); |
304 | } | |
305 | //------------------------------------------- | |
306 | //---------- Recluster cells ---------------- | |
307 | //------------------------------------------- | |
308 | ||
309 | else{ | |
310 | //------------------------------------------------------------------------------------- | |
311 | //Transform CaloCells into Digits | |
312 | //------------------------------------------------------------------------------------- | |
6060ed91 | 313 | |
314 | //In case of MC, first loop on the clusters and fill MC label to array | |
315 | //..................................................................... | |
316 | ||
317 | // for(Int_t j = 0; j < 24*48*11; j++) { | |
318 | // if(fCellLabels[j]!=-1) printf("Not well initialized cell %d, label %d\n",j,fCellLabels[j]) ; | |
319 | // } | |
320 | ||
321 | for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++) | |
322 | { | |
323 | AliVCluster *clus = event->GetCaloCluster(i); | |
324 | if(clus->IsEMCAL()){ | |
325 | ||
326 | Int_t label = clus->GetLabel(); | |
327 | UShort_t * index = clus->GetCellsAbsId() ; | |
328 | for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){ | |
329 | fCellLabels[index[icell]]=label; | |
330 | //printf("1) absID %d, label %d\n",index[icell], fCellLabels[index[icell]]); | |
331 | } | |
332 | nClustersOrg++; | |
333 | } | |
334 | } | |
335 | ||
336 | // Create digits | |
337 | //................. | |
6fdd30c4 | 338 | Int_t idigit = 0; |
339 | Int_t id = -1; | |
340 | Float_t amp = -1; | |
341 | Float_t time = -1; | |
342 | ||
343 | Double_t cellAmplitude = 0; | |
344 | Double_t cellTime = 0; | |
345 | Short_t cellNumber = 0; | |
346 | ||
347 | TTree *digitsTree = new TTree("digitstree","digitstree"); | |
348 | digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000); | |
349 | ||
6060ed91 | 350 | |
6fdd30c4 | 351 | for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) |
352 | { | |
353 | if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE) | |
354 | break; | |
355 | ||
356 | time = cellTime; | |
357 | amp = cellAmplitude; | |
358 | id = cellNumber; | |
359 | ||
5994e71f | 360 | Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; |
361 | fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta); | |
362 | fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); | |
363 | ||
364 | //Do not include bad channels found in analysis? | |
365 | if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && | |
366 | fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ | |
6060ed91 | 367 | fCellLabels[id]=-1; //reset the entry in the array for next event |
5994e71f | 368 | //printf("Remove channel %d\n",id); |
369 | continue; | |
370 | } | |
371 | ||
372 | //Recalibrate? | |
373 | if(fRecoUtils->IsRecalibrationOn()){ | |
374 | //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id); | |
375 | amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi); | |
376 | } | |
377 | ||
6060ed91 | 378 | //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary |
379 | new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); | |
380 | //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ; | |
381 | //else printf("2) Digit cell %d, no label, amp %f \n",id,amp) ; | |
382 | fCellLabels[id]=-1; //reset the entry in the array for next event | |
383 | ||
384 | //AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit); | |
385 | //digit->SetId(id); | |
386 | //digit->SetAmplitude(amp); | |
387 | //digit->SetTime(time); | |
388 | //digit->SetTimeR(time); | |
389 | //digit->SetIndexInList(idigit); | |
390 | //digit->SetType(AliEMCALDigit::kHG); | |
5994e71f | 391 | |
392 | //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit); | |
6fdd30c4 | 393 | idigit++; |
394 | } | |
395 | ||
396 | //Fill the tree with digits | |
397 | digitsTree->Fill(); | |
398 | ||
399 | //------------------------------------------------------------------------------------- | |
400 | //Do the clusterization | |
401 | //------------------------------------------------------------------------------------- | |
402 | TTree *clustersTree = new TTree("clustertree","clustertree"); | |
3fa4fb85 | 403 | |
6fdd30c4 | 404 | fClusterizer->SetInput(digitsTree); |
405 | fClusterizer->SetOutput(clustersTree); | |
406 | fClusterizer->Digits2Clusters(""); | |
6060ed91 | 407 | |
6fdd30c4 | 408 | //------------------------------------------------------------------------------------- |
409 | //Transform the recpoints into AliVClusters | |
410 | //------------------------------------------------------------------------------------- | |
411 | ||
412 | clustersTree->SetBranchStatus("*",0); //disable all branches | |
413 | clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need | |
414 | ||
415 | TBranch *branch = clustersTree->GetBranch("EMCALECARP"); | |
416 | branch->SetAddress(&fClusterArr); | |
417 | branch->GetEntry(0); | |
6060ed91 | 418 | |
6fdd30c4 | 419 | RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr); |
420 | ||
421 | //---CLEAN UP----- | |
422 | fClusterizer->Clear(); | |
423 | fDigitsArr ->Clear("C"); | |
424 | fClusterArr ->Delete(); // Do not Clear(), it leaks, why? | |
425 | ||
426 | clustersTree->Delete("all"); | |
427 | digitsTree ->Delete("all"); | |
6fdd30c4 | 428 | } |
429 | ||
385b7abf | 430 | //Recalculate track-matching for the new clusters, only with ESDs |
431 | if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr); | |
432 | ||
433 | ||
6fdd30c4 | 434 | //------------------------------------------------------------------------------------- |
435 | //Put the new clusters in the AOD list | |
436 | //------------------------------------------------------------------------------------- | |
437 | ||
438 | Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntries(); | |
6060ed91 | 439 | //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg); |
6fdd30c4 | 440 | for(Int_t i = 0; i < kNumberOfCaloClusters; i++){ |
441 | AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i); | |
3fa4fb85 | 442 | //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E()); |
385b7abf | 443 | |
444 | //Add matched track, if any, only with ESDs | |
445 | if(esdevent){ | |
446 | Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i); | |
447 | if(trackIndex >= 0){ | |
448 | newCluster->AddTrackMatched(event->GetTrack(trackIndex)); | |
449 | if(DebugLevel() > 1) | |
450 | Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i); | |
451 | } | |
452 | } | |
0821ba9f | 453 | |
5994e71f | 454 | //In case of new bad channels, recalculate distance to bad channels |
455 | if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){ | |
456 | //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel()); | |
457 | fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster); | |
458 | //printf("; after %f \n",newCluster->GetDistanceToBadChannel()); | |
459 | } | |
460 | ||
6060ed91 | 461 | // if(newCluster->GetNLabels()>0){ |
462 | // printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() ); | |
463 | // UShort_t * newindex = newCluster->GetCellsAbsId() ; | |
464 | // for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){ | |
465 | // printf("\t absID %d\n",newindex[icell]); | |
466 | // } | |
467 | // } | |
468 | // | |
469 | // printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E()); | |
470 | ||
0821ba9f | 471 | newCluster->SetID(i); |
6fdd30c4 | 472 | new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster); |
473 | } | |
474 | ||
475 | //---CLEAN UP----- | |
476 | fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why? | |
3fa4fb85 | 477 | } |
6fdd30c4 | 478 | |
479 | //_____________________________________________________________________ | |
bd9e8ebd | 480 | Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB() |
6fdd30c4 | 481 | { |
482 | //Access to OCDB stuff | |
bd9e8ebd | 483 | |
6fdd30c4 | 484 | AliVEvent * event = InputEvent(); |
bd9e8ebd | 485 | if (!event) |
6fdd30c4 | 486 | { |
385b7abf | 487 | Warning("AccessOCDB","Event not available!!!"); |
bd9e8ebd | 488 | return kFALSE; |
489 | } | |
490 | ||
491 | if (event->GetRunNumber()==fRun) | |
492 | return kTRUE; | |
493 | fRun = event->GetRunNumber(); | |
494 | ||
495 | if(DebugLevel() > 1 ) | |
496 | Info("AccessODCD()"," Begin"); | |
497 | ||
498 | fGeom = AliEMCALGeometry::GetInstance(fGeomName); | |
385b7abf | 499 | |
500 | ||
bd9e8ebd | 501 | AliCDBManager *cdb = AliCDBManager::Instance(); |
385b7abf | 502 | |
503 | ||
504 | if (fOCDBpath.Length()){ | |
6fdd30c4 | 505 | cdb->SetDefaultStorage(fOCDBpath.Data()); |
385b7abf | 506 | Info("AccessOCDB","Default storage %s",fOCDBpath.Data()); |
507 | } | |
508 | ||
bd9e8ebd | 509 | cdb->SetRun(event->GetRunNumber()); |
385b7abf | 510 | |
bd9e8ebd | 511 | // |
512 | // EMCAL from RAW OCDB | |
513 | if (fOCDBpath.Contains("alien:")) | |
514 | { | |
515 | cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB"); | |
516 | cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB"); | |
517 | } | |
385b7abf | 518 | |
bd9e8ebd | 519 | TString path = cdb->GetDefaultStorage()->GetBaseFolder(); |
385b7abf | 520 | |
bd9e8ebd | 521 | // init parameters: |
385b7abf | 522 | |
bd9e8ebd | 523 | //Get calibration parameters |
524 | if(!fCalibData) | |
525 | { | |
526 | AliCDBEntry *entry = (AliCDBEntry*) | |
527 | AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); | |
385b7abf | 528 | |
bd9e8ebd | 529 | if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); |
6fdd30c4 | 530 | } |
bd9e8ebd | 531 | |
532 | if(!fCalibData) | |
533 | AliFatal("Calibration parameters not found in CDB!"); | |
534 | ||
535 | //Get calibration parameters | |
536 | if(!fPedestalData) | |
6fdd30c4 | 537 | { |
bd9e8ebd | 538 | AliCDBEntry *entry = (AliCDBEntry*) |
539 | AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"); | |
385b7abf | 540 | |
bd9e8ebd | 541 | if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject(); |
6fdd30c4 | 542 | } |
bd9e8ebd | 543 | |
544 | if(!fPedestalData) | |
545 | AliFatal("Dead map not found in CDB!"); | |
6fdd30c4 | 546 | |
bd9e8ebd | 547 | InitClusterization(); |
385b7abf | 548 | |
6fdd30c4 | 549 | return kTRUE; |
6fdd30c4 | 550 | } |
551 | ||
552 | //________________________________________________________________________________________ | |
553 | void AliAnalysisTaskEMCALClusterize::InitClusterization() | |
554 | { | |
555 | //Select clusterization/unfolding algorithm and set all the needed parameters | |
556 | ||
bd9e8ebd | 557 | if (fJustUnfold){ |
558 | // init the unfolding afterburner | |
559 | delete fUnfolder; | |
560 | fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut()); | |
561 | return; | |
562 | } | |
563 | ||
564 | //First init the clusterizer | |
565 | delete fClusterizer; | |
566 | if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1) | |
567 | fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData); | |
568 | else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) | |
569 | fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData); | |
570 | else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) { | |
571 | AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData); | |
572 | clusterizer->SetNRowDiff(2); | |
573 | clusterizer->SetNColDiff(2); | |
574 | fClusterizer = clusterizer; | |
575 | } else{ | |
576 | AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag())); | |
577 | } | |
6fdd30c4 | 578 | |
bd9e8ebd | 579 | //Now set the parameters |
580 | fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() ); | |
581 | fClusterizer->SetECALogWeight ( fRecParam->GetW0() ); | |
582 | fClusterizer->SetMinECut ( fRecParam->GetMinECut() ); | |
583 | fClusterizer->SetUnfolding ( fRecParam->GetUnfold() ); | |
584 | fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() ); | |
585 | fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() ); | |
586 | fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() ); | |
587 | fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() ); | |
588 | fClusterizer->SetInputCalibrated ( kTRUE ); | |
2af35150 | 589 | fClusterizer->SetJustClusters ( kTRUE ); |
bd9e8ebd | 590 | //In case of unfolding after clusterization is requested, set the corresponding parameters |
591 | if(fRecParam->GetUnfold()){ | |
592 | Int_t i=0; | |
593 | for (i = 0; i < 8; i++) { | |
594 | fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i)); | |
595 | }//end of loop over parameters | |
596 | for (i = 0; i < 3; i++) { | |
597 | fClusterizer->SetPar5 (i, fRecParam->GetPar5(i)); | |
598 | fClusterizer->SetPar6 (i, fRecParam->GetPar6(i)); | |
599 | }//end of loop over parameters | |
6fdd30c4 | 600 | |
bd9e8ebd | 601 | fClusterizer->InitClusterUnfolding(); |
602 | }// to unfold | |
6fdd30c4 | 603 | } |
3fa4fb85 | 604 | |
6fdd30c4 | 605 | //________________________________________________________________________________________ |
606 | void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray) | |
607 | { | |
608 | // Restore clusters from recPoints | |
609 | // Cluster energy, global position, cells and their amplitude fractions are restored | |
610 | ||
611 | for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++) | |
612 | { | |
613 | AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i); | |
614 | ||
615 | const Int_t ncells = recPoint->GetMultiplicity(); | |
616 | Int_t ncells_true = 0; | |
617 | ||
618 | // cells and their amplitude fractions | |
619 | UShort_t absIds[ncells]; | |
620 | Double32_t ratios[ncells]; | |
621 | ||
622 | for (Int_t c = 0; c < ncells; c++) { | |
623 | AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]); | |
624 | ||
625 | absIds[ncells_true] = digit->GetId(); | |
626 | ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude(); | |
627 | ||
628 | if (ratios[ncells_true] > 0.001) ncells_true++; | |
629 | } | |
630 | ||
631 | if (ncells_true < 1) { | |
632 | Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells"); | |
633 | continue; | |
634 | } | |
635 | ||
636 | TVector3 gpos; | |
637 | Float_t g[3]; | |
638 | ||
639 | // calculate new cluster position | |
640 | recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr); | |
641 | recPoint->GetGlobalPosition(gpos); | |
642 | gpos.GetXYZ(g); | |
643 | ||
644 | // create a new cluster | |
645 | AliAODCaloCluster *clus = new AliAODCaloCluster(); | |
646 | clus->SetType(AliVCluster::kEMCALClusterv1); | |
647 | clus->SetE(recPoint->GetEnergy()); | |
648 | clus->SetPosition(g); | |
649 | clus->SetNCells(ncells_true); | |
650 | clus->SetCellsAbsId(absIds); | |
651 | clus->SetCellsAmplitudeFraction(ratios); | |
652 | clus->SetDispersion(recPoint->GetDispersion()); | |
653 | clus->SetChi2(-1); //not yet implemented | |
bd9e8ebd | 654 | clus->SetTOF(recPoint->GetTime()) ; //time-of-flight |
6fdd30c4 | 655 | clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima |
656 | Float_t elipAxis[2]; | |
657 | recPoint->GetElipsAxis(elipAxis); | |
658 | clus->SetM02(elipAxis[0]*elipAxis[0]) ; | |
659 | clus->SetM20(elipAxis[1]*elipAxis[1]) ; | |
5994e71f | 660 | clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); |
6060ed91 | 661 | |
662 | //MC | |
663 | Int_t parentMult = 0; | |
664 | Int_t *parentList = recPoint->GetParents(parentMult); | |
665 | clus->SetLabel(parentList, parentMult); | |
666 | ||
6fdd30c4 | 667 | clusArray->Add(clus); |
6fdd30c4 | 668 | } // recPoints loop |
6fdd30c4 | 669 | } |