]>
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 | **************************************************************************/ | |
6fdd30c4 | 15 | |
16 | //_________________________________________________________________________ | |
17 | // This analysis provides a new list of clusters to be used in other analysis | |
18 | // | |
19 | // Author: Gustavo Conesa Balbastre, | |
20 | // Adapted from analysis class from Deepa Thomas | |
21 | // | |
22 | // | |
23 | //_________________________________________________________________________ | |
24 | ||
25 | // --- Root --- | |
26 | #include "TString.h" | |
27 | #include "TRefArray.h" | |
28 | #include "TClonesArray.h" | |
29 | #include "TTree.h" | |
30 | #include "TGeoManager.h" | |
5994e71f | 31 | #include "TROOT.h" |
32 | #include "TInterpreter.h" | |
1f77507b | 33 | #include "TFile.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" |
3769e0cb | 48 | #include "AliOADBContainer.h" |
6fdd30c4 | 49 | |
50 | // --- EMCAL | |
6fdd30c4 | 51 | #include "AliEMCALAfterBurnerUF.h" |
52 | #include "AliEMCALGeometry.h" | |
53 | #include "AliEMCALClusterizerNxN.h" | |
54 | #include "AliEMCALClusterizerv1.h" | |
cfcbe5d2 | 55 | #include "AliEMCALClusterizerv2.h" |
6fdd30c4 | 56 | #include "AliEMCALRecPoint.h" |
57 | #include "AliEMCALDigit.h" | |
58 | #include "AliCaloCalibPedestal.h" | |
59 | #include "AliEMCALCalibData.h" | |
60 | ||
61 | #include "AliAnalysisTaskEMCALClusterize.h" | |
62 | ||
63 | ClassImp(AliAnalysisTaskEMCALClusterize) | |
64 | ||
6544055e | 65 | //______________________________________________________________________________ |
6fdd30c4 | 66 | AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) |
6544055e | 67 | : AliAnalysisTaskSE(name) |
68 | , fEvent(0) | |
3769e0cb | 69 | , fGeom(0), fGeomName("") |
6544055e | 70 | , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE) |
71 | , fCalibData(0), fPedestalData(0) | |
3769e0cb | 72 | , fOCDBpath(""), fAccessOCDB(kFALSE) |
6544055e | 73 | , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0) |
74 | , fRecParam(0), fClusterizer(0) | |
75 | , fUnfolder(0), fJustUnfold(kFALSE) | |
3769e0cb | 76 | , fOutputAODBranch(0), fOutputAODBranchName("") |
6544055e | 77 | , fFillAODFile(kTRUE), fFillAODHeader(0) |
78 | , fFillAODCaloCells(0), fRun(-1) | |
79 | , fRecoUtils(0), fConfigName("") | |
80 | , fCellLabels(), fCellSecondLabels(), fCellTime() | |
adad4ea9 | 81 | , fCellMatchdEta(), fCellMatchdPhi() |
3769e0cb | 82 | , fMaxEvent(0), fDoTrackMatching(kFALSE) |
83 | , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0) | |
84 | , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE) | |
85 | , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("") | |
86 | , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("") | |
6a8be5c3 | 87 | { |
6fdd30c4 | 88 | //ctor |
e3990982 | 89 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; |
90 | for(Int_t j = 0; j < 24*48*11; j++) | |
91 | { | |
6a8be5c3 | 92 | fCellLabels[j] = -1; |
93 | fCellSecondLabels[j] = -1; | |
adad4ea9 | 94 | fCellTime[j] = 0.; |
95 | fCellMatchdEta[j] = -999; | |
96 | fCellMatchdPhi[j] = -999; | |
6a8be5c3 | 97 | } |
a39f5b70 | 98 | |
a39f5b70 | 99 | |
6fdd30c4 | 100 | } |
101 | ||
6544055e | 102 | //______________________________________________________________ |
6fdd30c4 | 103 | AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() |
6544055e | 104 | : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize") |
105 | , fEvent(0) | |
3769e0cb | 106 | , fGeom(0), fGeomName("") |
6544055e | 107 | , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE) |
108 | , fCalibData(0), fPedestalData(0) | |
3769e0cb | 109 | , fOCDBpath(""), fAccessOCDB(kFALSE) |
6544055e | 110 | , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0) |
111 | , fRecParam(0), fClusterizer(0) | |
112 | , fUnfolder(0), fJustUnfold(kFALSE) | |
3769e0cb | 113 | , fOutputAODBranch(0), fOutputAODBranchName("") |
6544055e | 114 | , fFillAODFile(kTRUE), fFillAODHeader(0) |
115 | , fFillAODCaloCells(0), fRun(-1) | |
116 | , fRecoUtils(0), fConfigName("") | |
117 | , fCellLabels(), fCellSecondLabels(), fCellTime() | |
adad4ea9 | 118 | , fCellMatchdEta(), fCellMatchdPhi() |
3769e0cb | 119 | , fMaxEvent(0), fDoTrackMatching(kFALSE) |
120 | , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0) | |
121 | , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE) | |
9ac44255 | 122 | , fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("") |
3769e0cb | 123 | , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("") |
124 | ||
6fdd30c4 | 125 | { |
126 | // Constructor | |
e3990982 | 127 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; |
128 | for(Int_t j = 0; j < 24*48*11; j++) | |
129 | { | |
37d2296c | 130 | fCellLabels[j] = -1; |
131 | fCellSecondLabels[j] = -1; | |
adad4ea9 | 132 | fCellTime[j] = 0.; |
133 | fCellMatchdEta[j] = -999; | |
134 | fCellMatchdPhi[j] = -999; | |
37d2296c | 135 | } |
3769e0cb | 136 | |
6fdd30c4 | 137 | } |
138 | ||
5994e71f | 139 | |
6544055e | 140 | //_______________________________________________________________ |
6fdd30c4 | 141 | AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize() |
142 | { | |
143 | //dtor | |
144 | ||
145 | if (fDigitsArr){ | |
146 | fDigitsArr->Clear("C"); | |
542e1379 | 147 | delete fDigitsArr; |
6fdd30c4 | 148 | } |
149 | ||
150 | if (fClusterArr){ | |
151 | fClusterArr->Delete(); | |
542e1379 | 152 | delete fClusterArr; |
6fdd30c4 | 153 | } |
154 | ||
155 | if (fCaloClusterArr){ | |
156 | fCaloClusterArr->Delete(); | |
542e1379 | 157 | delete fCaloClusterArr; |
6fdd30c4 | 158 | } |
6a8be5c3 | 159 | |
385b7abf | 160 | if(fClusterizer) delete fClusterizer; |
161 | if(fUnfolder) delete fUnfolder; | |
162 | if(fRecoUtils) delete fRecoUtils; | |
6a8be5c3 | 163 | |
6fdd30c4 | 164 | } |
165 | ||
3769e0cb | 166 | //_______________________________________________ |
167 | void AliAnalysisTaskEMCALClusterize::AccessOADB() | |
168 | { | |
169 | // Set the AODB calibration, bad channels etc. parameters at least once | |
170 | // alignment matrices from OADB done in SetGeometryMatrices | |
171 | ||
172 | //Set it only once | |
173 | if(fOADBSet) return ; | |
174 | ||
175 | Int_t runnumber = InputEvent()->GetRunNumber() ; | |
176 | TString pass = GetPass(); | |
177 | ||
178 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data()); | |
179 | ||
180 | Int_t nSM = fGeom->GetNumberOfSuperModules(); | |
181 | ||
182 | // Bad map | |
183 | if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()) | |
184 | { | |
185 | AliOADBContainer *contBC=new AliOADBContainer(""); | |
186 | contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels"); | |
187 | ||
188 | TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber); | |
189 | ||
190 | if(arrayBC) | |
191 | { | |
192 | TObjArray * arrayBCpass = 0x0; | |
193 | if(pass!="")arrayBCpass = (TObjArray*)arrayBC->FindObject(pass); | |
194 | // There are no passes for simulation, in order to get the bad map put pass1 | |
195 | else arrayBCpass = (TObjArray*)arrayBC->FindObject("pass1"); | |
196 | ||
197 | if(arrayBCpass) | |
198 | { | |
199 | fRecoUtils->SwitchOnDistToBadChannelRecalculation(); | |
200 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n"); | |
201 | ||
202 | for (Int_t i=0; i<nSM; ++i) | |
203 | { | |
204 | TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i); | |
205 | ||
206 | if (hbm) | |
207 | delete hbm; | |
208 | ||
209 | hbm=(TH2I*)arrayBCpass->FindObject(Form("EMCALBadChannelMap_Mod%d",i)); | |
210 | ||
211 | if (!hbm) | |
212 | { | |
213 | AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i)); | |
214 | continue; | |
215 | } | |
216 | ||
217 | hbm->SetDirectory(0); | |
218 | fRecoUtils->SetEMCALChannelStatusMap(i,hbm); | |
219 | ||
220 | } // loop | |
221 | } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels 1\n"); // pass array | |
222 | } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels 2\n"); // run array | |
223 | } // Remove bad | |
224 | ||
225 | // Energy Recalibration | |
226 | if(fRecoUtils->IsRecalibrationOn()) | |
227 | { | |
228 | AliOADBContainer *contRF=new AliOADBContainer(""); | |
229 | ||
230 | contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib"); | |
231 | ||
232 | TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber); | |
233 | ||
234 | if(recal) | |
235 | { | |
236 | TObjArray *recalpass=(TObjArray*)recal->FindObject(pass); | |
237 | ||
238 | if(recalpass) | |
239 | { | |
240 | TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib"); | |
241 | ||
242 | if(recalib) | |
243 | { | |
244 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n"); | |
245 | for (Int_t i=0; i<nSM; ++i) | |
246 | { | |
247 | TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i); | |
248 | ||
249 | if (h) | |
250 | delete h; | |
251 | ||
252 | h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i)); | |
253 | ||
254 | if (!h) | |
255 | { | |
256 | AliError(Form("Could not load EMCALRecalFactors_SM%d",i)); | |
257 | continue; | |
258 | } | |
259 | ||
260 | h->SetDirectory(0); | |
261 | ||
262 | fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h); | |
263 | } // SM loop | |
264 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL 1\n"); // array ok | |
265 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL 2\n"); // array pass ok | |
266 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL 3\n"); // run number array ok | |
267 | ||
268 | // once set, apply run dependent corrections if requested | |
269 | fRecoUtils->SetRunDependentCorrections(runnumber); | |
270 | ||
271 | } // Recalibration on | |
272 | ||
273 | ||
274 | // Time Recalibration | |
275 | if(fRecoUtils->IsTimeRecalibrationOn()) | |
276 | { | |
277 | AliOADBContainer *contTRF=new AliOADBContainer(""); | |
278 | ||
f46af216 | 279 | contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib"); |
3769e0cb | 280 | |
281 | TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber); | |
282 | ||
283 | if(trecal) | |
284 | { | |
285 | TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass); | |
286 | ||
287 | if(trecalpass) | |
288 | { | |
f46af216 | 289 | printf("AliCalorimeterUtils::SetOADBParameters() - Time Recalibrate EMCAL \n"); |
290 | for (Int_t ibc = 0; ibc < 4; ++ibc) | |
3769e0cb | 291 | { |
f46af216 | 292 | TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc); |
293 | ||
294 | if (h) | |
295 | delete h; | |
296 | ||
297 | h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc)); | |
298 | ||
299 | if (!h) | |
3769e0cb | 300 | { |
f46af216 | 301 | AliError(Form("Could not load hAllTimeAvBC%d",ibc)); |
302 | continue; | |
303 | } | |
304 | ||
305 | h->SetDirectory(0); | |
306 | ||
307 | fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h); | |
308 | } // bunch crossing loop | |
309 | }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL 1\n"); // array pass ok | |
310 | }else printf("AliCalorimeterUtils::SetOADBParameters() - Do NOT recalibrate time EMCAL 2\n"); // run number array ok | |
3769e0cb | 311 | |
312 | } // Time recalibration on | |
313 | ||
314 | // Parameters already set once, so do not it again | |
315 | fOADBSet = kTRUE; | |
316 | ||
317 | } | |
318 | ||
6544055e | 319 | //_________________________________________________ |
320 | Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB() | |
321 | { | |
322 | //Access to OCDB stuff | |
1f77507b | 323 | |
6544055e | 324 | fEvent = InputEvent(); |
325 | if (!fEvent) | |
326 | { | |
327 | Warning("AccessOCDB","Event not available!!!"); | |
328 | return kFALSE; | |
1f77507b | 329 | } |
1f77507b | 330 | |
6544055e | 331 | if (fEvent->GetRunNumber()==fRun) |
332 | return kTRUE; | |
333 | fRun = fEvent->GetRunNumber(); | |
1f77507b | 334 | |
6544055e | 335 | if(DebugLevel() > 1 ) |
336 | printf("AliAnalysisTaksEMCALClusterize::AccessODCD() - Begin"); | |
3769e0cb | 337 | |
6544055e | 338 | AliCDBManager *cdb = AliCDBManager::Instance(); |
1f77507b | 339 | |
1f77507b | 340 | |
6544055e | 341 | if (fOCDBpath.Length()){ |
342 | cdb->SetDefaultStorage(fOCDBpath.Data()); | |
343 | printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data()); | |
344 | } | |
345 | ||
346 | cdb->SetRun(fEvent->GetRunNumber()); | |
1f77507b | 347 | |
348 | // | |
6544055e | 349 | // EMCAL from RAW OCDB |
350 | if (fOCDBpath.Contains("alien:")) | |
351 | { | |
352 | cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB"); | |
353 | cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB"); | |
354 | } | |
1f77507b | 355 | |
6544055e | 356 | TString path = cdb->GetDefaultStorage()->GetBaseFolder(); |
1f77507b | 357 | |
6544055e | 358 | // init parameters: |
1f77507b | 359 | |
6544055e | 360 | //Get calibration parameters |
361 | if(!fCalibData) | |
362 | { | |
363 | AliCDBEntry *entry = (AliCDBEntry*) | |
364 | AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); | |
365 | ||
366 | if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); | |
1f77507b | 367 | } |
6544055e | 368 | |
369 | if(!fCalibData) | |
370 | AliFatal("Calibration parameters not found in CDB!"); | |
371 | ||
372 | //Get calibration parameters | |
373 | if(!fPedestalData) | |
374 | { | |
375 | AliCDBEntry *entry = (AliCDBEntry*) | |
376 | AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"); | |
377 | ||
378 | if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject(); | |
1f77507b | 379 | } |
380 | ||
6544055e | 381 | if(!fPedestalData) |
382 | AliFatal("Dead map not found in CDB!"); | |
1f77507b | 383 | |
6544055e | 384 | return kTRUE; |
6fdd30c4 | 385 | } |
386 | ||
9ac44255 | 387 | //_____________________________________________________ |
6544055e | 388 | void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent() |
6fdd30c4 | 389 | { |
6544055e | 390 | // Get the input event, it can depend in embedded events what you want to get |
391 | // Also check if the quality of the event is good if not reject it | |
392 | ||
393 | fEvent = 0x0; | |
6a8be5c3 | 394 | |
220447cd | 395 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); |
396 | Int_t eventN = Entry(); | |
397 | if(aodIH) eventN = aodIH->GetReadEntry(); | |
398 | ||
6544055e | 399 | if (eventN > fMaxEvent) |
400 | return ; | |
f5036bcb | 401 | |
c50e6949 | 402 | //printf("Clusterizer --- Event %d-- \n",eventN); |
6544055e | 403 | |
6060ed91 | 404 | //Check if input event are embedded events |
405 | //If so, take output event | |
3769e0cb | 406 | if (aodIH && aodIH->GetMergeEvents()) |
407 | { | |
6544055e | 408 | fEvent = AODEvent(); |
1f77507b | 409 | |
410 | if(!aodIH->GetMergeEMCALCells()) | |
411 | AliFatal("Events merged but not EMCAL cells, check analysis settings!"); | |
412 | ||
3769e0cb | 413 | if(DebugLevel() > 1) |
414 | { | |
415 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n"); | |
416 | ||
417 | printf("\t InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(), | |
418 | InputEvent()->GetEMCALCells()->GetNumberOfCells()); | |
419 | ||
420 | printf("\t MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), | |
421 | aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()); | |
422 | ||
423 | for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) | |
424 | { | |
425 | AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl); | |
426 | if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E %f\n",icl,sigCluster->E()); | |
427 | } | |
428 | ||
429 | printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(), | |
430 | AODEvent()->GetEMCALCells()->GetNumberOfCells()); | |
6544055e | 431 | } |
6060ed91 | 432 | } |
3769e0cb | 433 | else |
434 | { | |
6544055e | 435 | fEvent = InputEvent(); |
1f77507b | 436 | if(fFillAODCaloCells) FillAODCaloCells(); |
437 | if(fFillAODHeader) FillAODHeader(); | |
6060ed91 | 438 | } |
439 | ||
3769e0cb | 440 | if (!fEvent) |
441 | { | |
3fa4fb85 | 442 | Error("UserExec","Event not available"); |
6544055e | 443 | return ; |
6fdd30c4 | 444 | } |
6a8be5c3 | 445 | |
cd2e4ce6 | 446 | //------------------------------------------------------------------------------------- |
6544055e | 447 | // Reject events if LED was firing, use only for LHC11a data |
448 | // Reject event if triggered by exotic cell and remove exotic cells if not triggered | |
cd2e4ce6 | 449 | //------------------------------------------------------------------------------------- |
6544055e | 450 | |
3769e0cb | 451 | if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; } |
6544055e | 452 | |
3769e0cb | 453 | if( IsExoticEvent() ) { fEvent = 0x0 ; return ; } |
cd2e4ce6 | 454 | |
f5036bcb | 455 | //Magic line to write events to AOD filem put after event rejection |
456 | AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile); | |
6544055e | 457 | |
458 | } | |
459 | ||
460 | //____________________________________________________ | |
461 | void AliAnalysisTaskEMCALClusterize::ClusterizeCells() | |
462 | { | |
463 | // Recluster calocells, transform them into digits, | |
464 | // feed the clusterizer with them and get new list of clusters | |
465 | ||
466 | //In case of MC, first loop on the clusters and fill MC label to array | |
467 | Int_t nClusters = fEvent->GetNumberOfCaloClusters(); | |
468 | Int_t nClustersOrg = 0; | |
f5036bcb | 469 | |
6544055e | 470 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); |
471 | if(aodIH && aodIH->GetEventToMerge()) //Embedding | |
472 | nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal | |
f5036bcb | 473 | |
6544055e | 474 | for (Int_t i = 0; i < nClusters; i++) |
475 | { | |
476 | AliVCluster *clus = 0; | |
477 | if(aodIH && aodIH->GetEventToMerge()) //Embedding | |
478 | clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal | |
479 | else | |
480 | clus = fEvent->GetCaloCluster(i); | |
5994e71f | 481 | |
6544055e | 482 | if(!clus) return; |
6a8be5c3 | 483 | |
f46af216 | 484 | if(clus->IsEMCAL()) |
485 | { | |
6544055e | 486 | Int_t label = clus->GetLabel(); |
487 | Int_t label2 = -1 ; | |
73a4ed8f | 488 | //printf("Org cluster E %f, Time %e, Id = ", clus->E(), clus->GetTOF() ); |
6544055e | 489 | if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ; |
490 | UShort_t * index = clus->GetCellsAbsId() ; | |
f46af216 | 491 | for(Int_t icell=0; icell < clus->GetNCells(); icell++ ) |
492 | { | |
6544055e | 493 | fCellLabels[index[icell]] = label; |
494 | fCellSecondLabels[index[icell]] = label2; | |
73a4ed8f | 495 | fCellTime[index[icell]] = clus->GetTOF(); |
496 | fCellMatchdEta[index[icell]] = clus->GetTrackDz(); | |
497 | fCellMatchdPhi[index[icell]] = clus->GetTrackDx(); | |
498 | //printf(" %d,", index[icell] ); | |
6fdd30c4 | 499 | } |
6544055e | 500 | nClustersOrg++; |
6fdd30c4 | 501 | } |
73a4ed8f | 502 | // printf("\n"); |
6544055e | 503 | } |
504 | ||
505 | // Transform CaloCells into Digits | |
506 | ||
507 | Int_t idigit = 0; | |
508 | Int_t id = -1; | |
509 | Float_t amp = -1; | |
510 | Double_t time = -1; | |
511 | ||
512 | AliVCaloCells *cells = fEvent->GetEMCALCells(); | |
3769e0cb | 513 | |
514 | TFile* file = new TFile("tmpTree.root","recreate"); // Trick to avoid annoying messages FIXME | |
6544055e | 515 | TTree *digitsTree = new TTree("digitstree","digitstree"); |
516 | digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000); | |
3769e0cb | 517 | |
a7e5a381 | 518 | Int_t bc = InputEvent()->GetBunchCrossNumber(); |
6544055e | 519 | for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) |
520 | { | |
6fdd30c4 | 521 | |
6544055e | 522 | // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell |
523 | id = cells->GetCellNumber(icell); | |
a7e5a381 | 524 | Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells); |
525 | ||
526 | // Do not include cells with too low energy, nor exotic cell | |
527 | if(amp < fRecParam->GetMinECut() ) accept = kFALSE; | |
528 | ||
529 | // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs. | |
3769e0cb | 530 | if (time*1e9 < 1.) |
531 | { | |
a7e5a381 | 532 | time = fCellTime[id]; |
f46af216 | 533 | //printf("cell %d time org %f - ",id, time*1.e9); |
a7e5a381 | 534 | fRecoUtils->RecalibrateCellTime(id,bc,time); |
f46af216 | 535 | //printf("recal %f\n",time*1.e9); |
a7e5a381 | 536 | } |
537 | ||
3769e0cb | 538 | if( accept && fRecoUtils->IsExoticCell(id,cells,bc)) |
539 | { | |
a7e5a381 | 540 | accept = kFALSE; |
541 | } | |
542 | ||
f46af216 | 543 | if( !accept ) |
544 | { | |
6544055e | 545 | fCellLabels[id] =-1; //reset the entry in the array for next event |
546 | fCellSecondLabels[id]=-1; //reset the entry in the array for next event | |
adad4ea9 | 547 | fCellTime[id] = 0.; |
548 | fCellMatchdEta[id] =-999; | |
549 | fCellMatchdPhi[id] =-999; | |
f46af216 | 550 | if( DebugLevel() > 2 ) |
551 | printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n", | |
6544055e | 552 | id,icell, cells->GetNumberOfCells(), amp, time*1.e9); |
553 | continue; | |
554 | } | |
b65291bf | 555 | |
6544055e | 556 | //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary |
557 | new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1); | |
558 | ||
559 | fCellLabels[id] =-1; //reset the entry in the array for next event | |
6fdd30c4 | 560 | |
6544055e | 561 | idigit++; |
562 | } | |
6a8be5c3 | 563 | |
6544055e | 564 | //Fill the tree with digits |
565 | digitsTree->Fill(); | |
6fdd30c4 | 566 | |
6544055e | 567 | //------------------------------------------------------------------------------------- |
568 | //Do the clusterization | |
569 | //------------------------------------------------------------------------------------- | |
570 | TTree *clustersTree = new TTree("clustertree","clustertree"); | |
571 | ||
572 | fClusterizer->SetInput(digitsTree); | |
573 | fClusterizer->SetOutput(clustersTree); | |
574 | fClusterizer->Digits2Clusters(""); | |
575 | ||
576 | //------------------------------------------------------------------------------------- | |
577 | //Transform the recpoints into AliVClusters | |
578 | //------------------------------------------------------------------------------------- | |
579 | ||
580 | clustersTree->SetBranchStatus("*",0); //disable all branches | |
581 | clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need | |
582 | ||
583 | TBranch *branch = clustersTree->GetBranch("EMCALECARP"); | |
584 | branch->SetAddress(&fClusterArr); | |
3769e0cb | 585 | |
6544055e | 586 | branch->GetEntry(0); |
3769e0cb | 587 | |
6544055e | 588 | RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr); |
589 | ||
3769e0cb | 590 | if(!fCaloClusterArr) |
591 | { | |
6544055e | 592 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast()); |
593 | return; | |
594 | } | |
595 | ||
3769e0cb | 596 | if( DebugLevel() > 0 ) |
597 | { | |
6544055e | 598 | printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast()); |
3bfc4732 | 599 | |
6544055e | 600 | if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){ |
601 | printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n", | |
602 | fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(), | |
603 | fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries()); | |
604 | } | |
605 | } | |
606 | ||
607 | //Reset the array with second labels for this event | |
608 | memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels)); | |
609 | ||
610 | //---CLEAN UP----- | |
611 | fClusterizer->Clear(); | |
612 | fDigitsArr ->Clear("C"); | |
613 | fClusterArr ->Delete(); // Do not Clear(), it leaks, why? | |
614 | ||
3769e0cb | 615 | clustersTree->Delete("all"); |
616 | digitsTree ->Delete("all"); | |
6544055e | 617 | |
3769e0cb | 618 | // Trick to avoid annoying messages FIXME |
619 | file->Close(); | |
620 | delete file; | |
6544055e | 621 | } |
a39f5b70 | 622 | |
73a4ed8f | 623 | //_____________________________________________________ |
6544055e | 624 | void AliAnalysisTaskEMCALClusterize::ClusterUnfolding() |
625 | { | |
626 | // Take the event clusters and unfold them | |
627 | ||
628 | AliVCaloCells *cells = fEvent->GetEMCALCells(); | |
629 | Double_t cellAmplitude = 0; | |
630 | Double_t cellTime = 0; | |
631 | Short_t cellNumber = 0; | |
632 | Int_t nClustersOrg = 0; | |
633 | ||
634 | // Fill the array with the EMCAL clusters, copy them | |
635 | for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++) | |
636 | { | |
637 | AliVCluster *clus = fEvent->GetCaloCluster(i); | |
638 | if(clus->IsEMCAL()){ | |
5994e71f | 639 | |
6544055e | 640 | //recalibrate/remove bad channels/etc if requested |
641 | if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){ | |
5994e71f | 642 | continue; |
6544055e | 643 | } |
3bfc4732 | 644 | |
6544055e | 645 | if(fRecoUtils->IsRecalibrationOn()){ |
646 | ||
647 | //Calibrate cluster | |
648 | fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells); | |
649 | ||
650 | //CalibrateCells | |
651 | for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) | |
652 | { | |
653 | if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE) | |
654 | break; | |
655 | ||
656 | Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; | |
657 | fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta); | |
658 | fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); | |
659 | ||
660 | //Do not include bad channels found in analysis? | |
661 | if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && | |
662 | fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ | |
663 | fCellLabels[cellNumber] =-1; //reset the entry in the array for next event | |
664 | fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event | |
adad4ea9 | 665 | fCellTime[cellNumber] = 0.; |
666 | fCellMatchdEta[cellNumber] =-999; | |
667 | fCellMatchdPhi[cellNumber] =-999; | |
6544055e | 668 | continue; |
669 | } | |
670 | ||
671 | cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime); | |
672 | ||
673 | }// cells loop | |
674 | }// recalibrate | |
5994e71f | 675 | |
6544055e | 676 | //Cast to ESD or AOD, needed to create the cluster array |
677 | AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus); | |
678 | AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus); | |
679 | if (esdCluster){ | |
680 | fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) ); | |
681 | }//ESD | |
682 | else if(aodCluster){ | |
683 | fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) ); | |
684 | }//AOD | |
685 | else | |
686 | Warning("UserExec()"," - Wrong CaloCluster type?"); | |
687 | nClustersOrg++; | |
6a8be5c3 | 688 | } |
6fdd30c4 | 689 | } |
690 | ||
6544055e | 691 | //Do the unfolding |
692 | fUnfolder->UnfoldClusters(fCaloClusterArr, cells); | |
6a8be5c3 | 693 | |
6544055e | 694 | //CLEAN-UP |
695 | fUnfolder->Clear(); | |
385b7abf | 696 | |
6544055e | 697 | } |
698 | ||
699 | //_____________________________________________________ | |
700 | void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() | |
701 | { | |
702 | // Put calo cells in standard branch | |
703 | AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells()); | |
704 | Int_t nEMcell = eventEMcells.GetNumberOfCells() ; | |
6fdd30c4 | 705 | |
6544055e | 706 | AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells()); |
707 | aodEMcells.CreateContainer(nEMcell); | |
708 | aodEMcells.SetType(AliVCaloCells::kEMCALCell); | |
709 | Double_t calibFactor = 1.; | |
710 | for (Int_t iCell = 0; iCell < nEMcell; iCell++) { | |
711 | Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; | |
712 | fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); | |
713 | fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); | |
385b7abf | 714 | |
6544055e | 715 | if(fRecoUtils->IsRecalibrationOn()){ |
716 | calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi); | |
385b7abf | 717 | } |
0821ba9f | 718 | |
6544055e | 719 | if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad |
720 | aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor); | |
721 | } | |
722 | else { | |
723 | aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0); | |
5994e71f | 724 | } |
6fdd30c4 | 725 | } |
6544055e | 726 | aodEMcells.Sort(); |
6fdd30c4 | 727 | |
6544055e | 728 | } |
6fdd30c4 | 729 | |
6544055e | 730 | //__________________________________________________ |
731 | void AliAnalysisTaskEMCALClusterize::FillAODHeader() | |
6fdd30c4 | 732 | { |
6544055e | 733 | //Put event header information in standard AOD branch |
e9dd2d80 | 734 | |
6544055e | 735 | AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent); |
736 | AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent); | |
e9dd2d80 | 737 | |
6544055e | 738 | Double_t pos[3] ; |
739 | Double_t covVtx[6]; | |
740 | for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.; | |
741 | ||
742 | AliAODHeader* header = AODEvent()->GetHeader(); | |
743 | header->SetRunNumber(fEvent->GetRunNumber()); | |
744 | ||
745 | if(esdevent){ | |
746 | TTree* tree = fInputHandler->GetTree(); | |
747 | if (tree) { | |
748 | TFile* file = tree->GetCurrentFile(); | |
749 | if (file) header->SetESDFileName(file->GetName()); | |
750 | } | |
751 | } | |
752 | else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName()); | |
e9dd2d80 | 753 | |
6544055e | 754 | header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber()); |
755 | header->SetOrbitNumber(fEvent->GetOrbitNumber()); | |
756 | header->SetPeriodNumber(fEvent->GetPeriodNumber()); | |
757 | header->SetEventType(fEvent->GetEventType()); | |
385b7abf | 758 | |
6544055e | 759 | //Centrality |
760 | if(fEvent->GetCentrality()){ | |
761 | header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality()))); | |
762 | } | |
763 | else{ | |
764 | header->SetCentrality(0); | |
765 | } | |
385b7abf | 766 | |
6544055e | 767 | //Trigger |
768 | header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection | |
769 | if (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses()); | |
770 | else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses()); | |
771 | header->SetTriggerMask(fEvent->GetTriggerMask()); | |
772 | header->SetTriggerCluster(fEvent->GetTriggerCluster()); | |
773 | if(esdevent){ | |
774 | header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs()); | |
775 | header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs()); | |
776 | header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs()); | |
777 | } | |
778 | else if (aodevent){ | |
779 | header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs()); | |
780 | header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs()); | |
781 | header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs()); | |
782 | } | |
e9dd2d80 | 783 | |
6544055e | 784 | header->SetMagneticField(fEvent->GetMagneticField()); |
785 | //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); | |
e9dd2d80 | 786 | |
6544055e | 787 | header->SetZDCN1Energy(fEvent->GetZDCN1Energy()); |
788 | header->SetZDCP1Energy(fEvent->GetZDCP1Energy()); | |
789 | header->SetZDCN2Energy(fEvent->GetZDCN2Energy()); | |
790 | header->SetZDCP2Energy(fEvent->GetZDCP2Energy()); | |
791 | header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1)); | |
385b7abf | 792 | |
6544055e | 793 | Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()}; |
794 | Float_t diamcov[3]; | |
795 | fEvent->GetDiamondCovXY(diamcov); | |
796 | header->SetDiamond(diamxy,diamcov); | |
797 | if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ()); | |
798 | else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ()); | |
e9dd2d80 | 799 | |
bd9e8ebd | 800 | // |
6544055e | 801 | // |
802 | Int_t nVertices = 1 ;/* = prim. vtx*/; | |
803 | Int_t nCaloClus = fEvent->GetNumberOfCaloClusters(); | |
e9dd2d80 | 804 | |
6544055e | 805 | AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0); |
385b7abf | 806 | |
6544055e | 807 | // Access to the AOD container of vertices |
808 | TClonesArray &vertices = *(AODEvent()->GetVertices()); | |
809 | Int_t jVertices=0; | |
385b7abf | 810 | |
6544055e | 811 | // Add primary vertex. The primary tracks will be defined |
812 | // after the loops on the composite objects (V0, cascades, kinks) | |
813 | fEvent->GetPrimaryVertex()->GetXYZ(pos); | |
814 | Float_t chi = 0; | |
815 | if (esdevent){ | |
816 | esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx); | |
817 | chi = esdevent->GetPrimaryVertex()->GetChi2toNDF(); | |
818 | } | |
819 | else if (aodevent){ | |
820 | aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx); | |
821 | chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD? | |
6fdd30c4 | 822 | } |
bd9e8ebd | 823 | |
6544055e | 824 | AliAODVertex * primary = new(vertices[jVertices++]) |
825 | AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary); | |
826 | primary->SetName(fEvent->GetPrimaryVertex()->GetName()); | |
827 | primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle()); | |
e9dd2d80 | 828 | |
6544055e | 829 | } |
830 | ||
3769e0cb | 831 | //_______________________________________________ |
832 | TString AliAnalysisTaskEMCALClusterize::GetPass() | |
833 | { | |
834 | // Get passx from filename. | |
835 | ||
836 | if (!AliAnalysisManager::GetAnalysisManager()->GetTree()) | |
837 | { | |
838 | AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n"); | |
839 | return TString(""); | |
840 | } | |
841 | ||
842 | if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()) | |
843 | { | |
844 | AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n"); | |
845 | return TString(""); | |
846 | } | |
847 | ||
848 | TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); | |
849 | if (pass.Contains("ass1")) return TString("pass1"); | |
850 | else if (pass.Contains("ass2")) return TString("pass2"); | |
851 | else if (pass.Contains("ass3")) return TString("pass3"); | |
852 | ||
853 | // No condition fullfilled, give a default value | |
854 | printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n"); | |
855 | return TString(""); | |
856 | ||
857 | } | |
858 | ||
6544055e | 859 | //_________________________________________ |
860 | void AliAnalysisTaskEMCALClusterize::Init() | |
861 | { | |
862 | //Init analysis with configuration macro if available | |
863 | ||
3769e0cb | 864 | fOADBSet = kFALSE; |
865 | if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ; | |
866 | ||
867 | fDigitsArr = new TClonesArray("AliEMCALDigit",12000); | |
868 | fClusterArr = new TObjArray(10000); | |
869 | fCaloClusterArr = new TObjArray(10000); | |
870 | fBranchNames = "ESD:AliESDHeader.,EMCALCells."; | |
871 | ||
872 | if(!fRecParam) fRecParam = new AliEMCALRecParam; | |
873 | if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils(); | |
874 | ||
875 | if(fMaxEvent <= 0) fMaxEvent = 1000000000; | |
876 | if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005; | |
877 | if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001; | |
878 | ||
879 | if (fOCDBpath == "") fOCDBpath = "raw://" ; | |
880 | if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ; | |
881 | ||
882 | if(gROOT->LoadMacro(fConfigName) >=0) | |
883 | { | |
6544055e | 884 | printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data()); |
885 | AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()"); | |
886 | fGeomName = clus->fGeomName; | |
887 | fLoadGeomMatrices = clus->fLoadGeomMatrices; | |
888 | fOCDBpath = clus->fOCDBpath; | |
889 | fAccessOCDB = clus->fAccessOCDB; | |
890 | fRecParam = clus->fRecParam; | |
891 | fJustUnfold = clus->fJustUnfold; | |
892 | fFillAODFile = clus->fFillAODFile; | |
893 | fRecoUtils = clus->fRecoUtils; | |
894 | fConfigName = clus->fConfigName; | |
895 | fMaxEvent = clus->fMaxEvent; | |
896 | fDoTrackMatching = clus->fDoTrackMatching; | |
897 | fOutputAODBranchName = clus->fOutputAODBranchName; | |
e3990982 | 898 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ; |
6fdd30c4 | 899 | } |
e9dd2d80 | 900 | |
9ac44255 | 901 | // Init geometry, I do not like much to do it like this ... |
3769e0cb | 902 | if(fImportGeometryFromFile && !gGeoManager) |
903 | { | |
904 | if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/PWGGA/EMCALTasks/macros/geometry.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root" | |
905 | printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data()); | |
906 | TGeoManager::Import(fImportGeometryFilePath) ; | |
9ac44255 | 907 | } |
3769e0cb | 908 | |
6544055e | 909 | } |
6fdd30c4 | 910 | |
6544055e | 911 | //_______________________________________________________ |
6fdd30c4 | 912 | void AliAnalysisTaskEMCALClusterize::InitClusterization() |
913 | { | |
914 | //Select clusterization/unfolding algorithm and set all the needed parameters | |
915 | ||
bd9e8ebd | 916 | if (fJustUnfold){ |
917 | // init the unfolding afterburner | |
918 | delete fUnfolder; | |
1759f743 | 919 | fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut()); |
bd9e8ebd | 920 | return; |
6a8be5c3 | 921 | } |
922 | ||
bd9e8ebd | 923 | //First init the clusterizer |
924 | delete fClusterizer; | |
925 | if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1) | |
926 | fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData); | |
cfcbe5d2 | 927 | else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) |
928 | fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData); | |
212ac797 | 929 | else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){ |
bd9e8ebd | 930 | fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData); |
212ac797 | 931 | fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); |
932 | fClusterizer->SetNColDiff(fRecParam->GetNColDiff()); | |
d390e7cb | 933 | } else { |
bd9e8ebd | 934 | AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag())); |
935 | } | |
6a8be5c3 | 936 | |
bd9e8ebd | 937 | //Now set the parameters |
938 | fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() ); | |
939 | fClusterizer->SetECALogWeight ( fRecParam->GetW0() ); | |
940 | fClusterizer->SetMinECut ( fRecParam->GetMinECut() ); | |
941 | fClusterizer->SetUnfolding ( fRecParam->GetUnfold() ); | |
942 | fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() ); | |
943 | fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() ); | |
944 | fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() ); | |
945 | fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() ); | |
946 | fClusterizer->SetInputCalibrated ( kTRUE ); | |
2af35150 | 947 | fClusterizer->SetJustClusters ( kTRUE ); |
bd9e8ebd | 948 | //In case of unfolding after clusterization is requested, set the corresponding parameters |
949 | if(fRecParam->GetUnfold()){ | |
950 | Int_t i=0; | |
951 | for (i = 0; i < 8; i++) { | |
952 | fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i)); | |
953 | }//end of loop over parameters | |
954 | for (i = 0; i < 3; i++) { | |
955 | fClusterizer->SetPar5 (i, fRecParam->GetPar5(i)); | |
956 | fClusterizer->SetPar6 (i, fRecParam->GetPar6(i)); | |
957 | }//end of loop over parameters | |
6fdd30c4 | 958 | |
bd9e8ebd | 959 | fClusterizer->InitClusterUnfolding(); |
6544055e | 960 | |
bd9e8ebd | 961 | }// to unfold |
6fdd30c4 | 962 | } |
3fa4fb85 | 963 | |
6544055e | 964 | //_________________________________________________ |
965 | void AliAnalysisTaskEMCALClusterize::InitGeometry() | |
966 | { | |
3769e0cb | 967 | // Init geometry and set the geometry matrix, for the first event, skip the rest |
6544055e | 968 | // Also set once the run dependent calibrations |
969 | ||
3769e0cb | 970 | if(fGeomMatrixSet) return; |
971 | ||
972 | Int_t runnumber = InputEvent()->GetRunNumber() ; | |
973 | if (!fGeom) | |
974 | { | |
975 | if(fGeomName=="") | |
976 | { | |
977 | if (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1"; | |
978 | else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1"; | |
979 | else fGeomName = "EMCAL_COMPLETE12SMV1"; | |
980 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber); | |
981 | } | |
6544055e | 982 | |
3769e0cb | 983 | fGeom = AliEMCALGeometry::GetInstance(fGeomName); |
6544055e | 984 | |
3769e0cb | 985 | if(fDebug > 0) |
986 | { | |
987 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber); | |
988 | if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices"); | |
989 | printf("\n"); | |
990 | } | |
991 | } // geometry pointer did not exist before | |
992 | ||
993 | if(fLoadGeomMatrices) | |
994 | { | |
995 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n"); | |
996 | ||
997 | // OADB if available | |
998 | AliOADBContainer emcGeoMat("AliEMCALgeo"); | |
999 | emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo"); | |
1000 | TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices"); | |
1001 | ||
1002 | ||
1003 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++) | |
1004 | { | |
1005 | ||
1006 | if (!fGeomMatrix[mod]) // Get it from OADB | |
1007 | { | |
1008 | if(fDebug > 1 ) | |
1009 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n", | |
1010 | mod,((TGeoHMatrix*) matEMCAL->At(mod))); | |
1011 | //((TGeoHMatrix*) matEMCAL->At(mod))->Print(); | |
1012 | ||
1013 | fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ; | |
1014 | } | |
1015 | ||
1016 | if(fGeomMatrix[mod]) | |
1017 | { | |
1018 | if(DebugLevel() > 1) | |
1019 | fGeomMatrix[mod]->Print(); | |
1020 | ||
1021 | fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ; | |
1022 | } | |
1023 | ||
1024 | fGeomMatrixSet=kTRUE; | |
1025 | ||
1026 | }//SM loop | |
1027 | }//Load matrices | |
1028 | else if(!gGeoManager) | |
1029 | { | |
1030 | printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data"); | |
1031 | //Still not implemented in AOD, just a workaround to be able to work at least with ESDs | |
1032 | if(!strcmp(fEvent->GetName(),"AliAODEvent")) | |
1033 | { | |
1034 | if(DebugLevel() > 1) | |
1035 | Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs."); | |
1036 | }//AOD | |
1037 | else | |
1038 | { | |
1039 | if(DebugLevel() > 1) | |
1040 | printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices."); | |
1041 | ||
1042 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ; | |
1043 | ||
1044 | if(!esd) | |
1045 | { | |
1046 | Error("InitGeometry"," - This event does not contain ESDs?"); | |
1047 | AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE); | |
1048 | return; | |
1049 | } | |
1050 | ||
1051 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++) | |
1052 | { | |
1053 | if(DebugLevel() > 1) | |
1054 | esd->GetEMCALMatrix(mod)->Print(); | |
1055 | ||
1056 | if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ; | |
1057 | ||
1058 | } | |
1059 | ||
1060 | fGeomMatrixSet=kTRUE; | |
1061 | ||
1062 | }//ESD | |
1063 | }//Load matrices from Data | |
1064 | ||
6544055e | 1065 | |
1066 | } | |
1067 | ||
6544055e | 1068 | //____________________________________________________ |
1069 | Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent() | |
1070 | { | |
1071 | ||
1072 | // Check if event is exotic, get an exotic cell and compare with triggered patch | |
1073 | // If there is a match remove event ... to be completed, filled with something provisional | |
1074 | ||
1075 | if(!fRemoveExoticEvents) return kFALSE; | |
1076 | ||
1077 | // Loop on cells | |
1078 | AliVCaloCells * cells = fEvent->GetEMCALCells(); | |
1079 | Float_t totCellE = 0; | |
a7e5a381 | 1080 | Int_t bc = InputEvent()->GetBunchCrossNumber(); |
6544055e | 1081 | for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++){ |
1082 | ||
1083 | Float_t ecell = 0 ; | |
1084 | Double_t tcell = 0 ; | |
1085 | ||
1086 | Int_t absID = cells->GetCellNumber(icell); | |
a7e5a381 | 1087 | Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells); |
1088 | if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell; | |
6544055e | 1089 | } |
1090 | ||
1091 | // TString triggerclasses = ""; | |
1092 | // if(esdevent) triggerclasses = esdevent ->GetFiredTriggerClasses(); | |
1093 | // else triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses(); | |
1094 | // // | |
1095 | // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4); | |
1096 | // AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE); | |
1097 | // return; | |
1098 | // | |
1099 | ||
1100 | //printf("TotE cell %f\n",totCellE); | |
1101 | if(totCellE < 1) return kTRUE; | |
1102 | ||
1103 | return kFALSE; | |
1104 | ||
1105 | } | |
1106 | ||
3769e0cb | 1107 | //________________________________________________________________ |
1108 | Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run) | |
6544055e | 1109 | { |
1110 | //Check if event is LED | |
1111 | ||
1112 | if(!fRemoveLEDEvents) return kFALSE; | |
3769e0cb | 1113 | |
1114 | //check events of LHC11a period | |
1115 | if(run < 140000 || run > 146860) return kFALSE ; | |
1116 | ||
6544055e | 1117 | // Count number of cells with energy larger than 0.1 in SM3, cut on this number |
1118 | Int_t ncellsSM3 = 0; | |
1119 | Int_t ncellsSM4 = 0; | |
1120 | AliVCaloCells * cells = fEvent->GetEMCALCells(); | |
3769e0cb | 1121 | for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) |
1122 | { | |
6544055e | 1123 | if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++; |
1124 | if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==4) ncellsSM4++; | |
1125 | } | |
1126 | ||
1127 | TString triggerclasses = ""; | |
1128 | ||
1129 | AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent); | |
1130 | if(esdevent) triggerclasses = esdevent ->GetFiredTriggerClasses(); | |
1131 | else triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses(); | |
1132 | ||
1133 | Int_t ncellcut = 21; | |
1134 | if(triggerclasses.Contains("EMC")) ncellcut = 35; | |
1135 | ||
3769e0cb | 1136 | if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 ) |
1137 | { | |
6544055e | 1138 | printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4); |
1139 | AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE); | |
1140 | return kTRUE; | |
1141 | } | |
1142 | ||
1143 | return kFALSE; | |
1144 | ||
1145 | } | |
1146 | ||
1147 | //______________________________________________________________________________ | |
1148 | void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, | |
1149 | TObjArray *recPoints, | |
1150 | TObjArray *clusArray) | |
6fdd30c4 | 1151 | { |
1152 | // Restore clusters from recPoints | |
1153 | // Cluster energy, global position, cells and their amplitude fractions are restored | |
6a8be5c3 | 1154 | Int_t j = 0; |
6fdd30c4 | 1155 | for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++) |
1156 | { | |
1157 | AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i); | |
1158 | ||
1159 | const Int_t ncells = recPoint->GetMultiplicity(); | |
f6e45fe0 | 1160 | Int_t ncellsTrue = 0; |
6fdd30c4 | 1161 | |
6a8be5c3 | 1162 | if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue; |
1163 | ||
6fdd30c4 | 1164 | // cells and their amplitude fractions |
1165 | UShort_t absIds[ncells]; | |
1166 | Double32_t ratios[ncells]; | |
1167 | ||
f6e45fe0 | 1168 | //For later check embedding |
1169 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); | |
f053064e | 1170 | |
1171 | Float_t clusterE = 0; | |
f46af216 | 1172 | for (Int_t c = 0; c < ncells; c++) |
1173 | { | |
6fdd30c4 | 1174 | AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]); |
1175 | ||
f6e45fe0 | 1176 | absIds[ncellsTrue] = digit->GetId(); |
1177 | ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude(); | |
6544055e | 1178 | |
d390e7cb | 1179 | // In case of unfolding, remove digits with unfolded energy too low |
afaaef51 | 1180 | if(fSelectCell){ |
f46af216 | 1181 | if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac) |
1182 | { | |
1183 | if(DebugLevel() > 1) | |
1184 | { | |
6544055e | 1185 | printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n", |
d390e7cb | 1186 | recPoint->GetEnergiesList()[c],digit->GetAmplitude()); |
1187 | } | |
1188 | ||
1189 | continue; | |
1190 | ||
1191 | } // if cuts | |
1192 | }// Select cells | |
6544055e | 1193 | |
d390e7cb | 1194 | //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected |
1195 | ncellsTrue++; | |
1196 | clusterE +=recPoint->GetEnergiesList()[c]; | |
f053064e | 1197 | |
f6e45fe0 | 1198 | // In case of embedding, fill ratio with amount of signal, |
f46af216 | 1199 | if (aodIH && aodIH->GetMergeEvents()) |
1200 | { | |
f6e45fe0 | 1201 | //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells(); |
1202 | AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells(); | |
1203 | AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells(); | |
1204 | ||
1205 | Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]); | |
1206 | //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]); | |
1207 | Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]); | |
1208 | //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude()); | |
1209 | ||
1210 | if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude; | |
1211 | //printf("\t \t ratio %f\n",ratios[ncellsTrue]); | |
1212 | ||
1213 | }//Embedding | |
1214 | ||
f6e45fe0 | 1215 | }// cluster cell loop |
6fdd30c4 | 1216 | |
f46af216 | 1217 | if (ncellsTrue < 1) |
1218 | { | |
f053064e | 1219 | if (DebugLevel() > 1) |
1220 | printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n", | |
1221 | recPoint->GetEnergy(), ncells); | |
1222 | continue; | |
1223 | } | |
1224 | ||
1225 | //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue); | |
1226 | ||
f46af216 | 1227 | if(clusterE < fRecParam->GetClusteringThreshold()) |
1228 | { | |
f053064e | 1229 | if (DebugLevel()>1) |
1230 | printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE); | |
6fdd30c4 | 1231 | continue; |
1232 | } | |
1233 | ||
1234 | TVector3 gpos; | |
1235 | Float_t g[3]; | |
1236 | ||
1237 | // calculate new cluster position | |
1238 | recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr); | |
1239 | recPoint->GetGlobalPosition(gpos); | |
1240 | gpos.GetXYZ(g); | |
1241 | ||
1242 | // create a new cluster | |
6a8be5c3 | 1243 | (*clusArray)[j] = new AliAODCaloCluster() ; |
1244 | AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ; | |
1245 | j++; | |
6fdd30c4 | 1246 | clus->SetType(AliVCluster::kEMCALClusterv1); |
6544055e | 1247 | clus->SetE(clusterE); |
6fdd30c4 | 1248 | clus->SetPosition(g); |
f6e45fe0 | 1249 | clus->SetNCells(ncellsTrue); |
6fdd30c4 | 1250 | clus->SetCellsAbsId(absIds); |
1251 | clus->SetCellsAmplitudeFraction(ratios); | |
6fdd30c4 | 1252 | clus->SetChi2(-1); //not yet implemented |
bd9e8ebd | 1253 | clus->SetTOF(recPoint->GetTime()) ; //time-of-flight |
6fdd30c4 | 1254 | clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima |
5994e71f | 1255 | clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); |
b65291bf | 1256 | |
f46af216 | 1257 | if(ncells == ncellsTrue) |
1258 | { | |
f053064e | 1259 | Float_t elipAxis[2]; |
1260 | recPoint->GetElipsAxis(elipAxis); | |
1261 | clus->SetM02(elipAxis[0]*elipAxis[0]) ; | |
1262 | clus->SetM20(elipAxis[1]*elipAxis[1]) ; | |
1263 | clus->SetDispersion(recPoint->GetDispersion()); | |
1264 | } | |
f46af216 | 1265 | else if(fSelectCell) |
1266 | { | |
afaaef51 | 1267 | // In case some cells rejected, in unfolding case, recalculate |
1268 | // shower shape parameters and position | |
f46af216 | 1269 | //if(DebugLevel() > 1) |
1270 | printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue); | |
f053064e | 1271 | AliVCaloCells* cells = 0x0; |
1272 | if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells(); | |
1273 | else cells = InputEvent()->GetEMCALCells(); | |
1274 | fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus); | |
1275 | fRecoUtils->RecalculateClusterPID(clus); | |
afaaef51 | 1276 | fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus); |
f053064e | 1277 | |
f053064e | 1278 | } |
6060ed91 | 1279 | |
1280 | //MC | |
b65291bf | 1281 | Int_t parentMult = 0; |
6060ed91 | 1282 | Int_t *parentList = recPoint->GetParents(parentMult); |
1283 | clus->SetLabel(parentList, parentMult); | |
1284 | ||
37d2296c | 1285 | //Write the second major contributor to each MC cluster. |
1286 | Int_t iNewLabel ; | |
f46af216 | 1287 | for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ) |
1288 | { | |
37d2296c | 1289 | |
1290 | Int_t idCell = clus->GetCellAbsId(iLoopCell) ; | |
6ffa213f | 1291 | if(idCell>=0){ |
1292 | iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label. | |
1293 | for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ ) | |
1294 | { | |
1295 | if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label. | |
1296 | if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0; | |
1297 | } | |
1298 | if (iNewLabel == 1) | |
1299 | { | |
937b6c91 | 1300 | Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ; |
f46af216 | 1301 | for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ) |
1302 | { | |
6ffa213f | 1303 | newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ; |
1304 | } | |
6a8be5c3 | 1305 | |
6ffa213f | 1306 | newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ; |
1307 | clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ; | |
abf0860a | 1308 | delete [] newLabelArray; |
6ffa213f | 1309 | } |
1310 | }//positive cell number | |
37d2296c | 1311 | } |
1312 | ||
6fdd30c4 | 1313 | } // recPoints loop |
6544055e | 1314 | |
1315 | } | |
1316 | ||
3769e0cb | 1317 | |
6544055e | 1318 | //____________________________________________________________ |
1319 | void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects() | |
1320 | { | |
1321 | // Init geometry, create list of output clusters | |
1322 | ||
3769e0cb | 1323 | if(fOutputAODBranchName.Length()!=0) |
1324 | { | |
6544055e | 1325 | fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0); |
1326 | fOutputAODBranch->SetName(fOutputAODBranchName); | |
b65291bf | 1327 | //fOutputAODBranch->SetOwner(kFALSE); |
6544055e | 1328 | AddAODBranch("TClonesArray", &fOutputAODBranch); |
1329 | } | |
3769e0cb | 1330 | else |
1331 | { | |
6544055e | 1332 | AliFatal("fOutputAODBranchName not set\n"); |
1333 | } | |
b65291bf | 1334 | |
1335 | //PostData(0,fOutputAODBranch); | |
1336 | ||
6fdd30c4 | 1337 | } |
6544055e | 1338 | |
1339 | //_______________________________________________________ | |
1340 | void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) | |
1341 | { | |
1342 | // Main loop | |
1343 | // Called for each event | |
3769e0cb | 1344 | |
6544055e | 1345 | //Remove the contents of output list set in the previous event |
1346 | fOutputAODBranch->Clear("C"); | |
1347 | ||
1348 | LoadBranches(); | |
1349 | ||
1350 | //Get the event, do some checks and settings | |
1351 | CheckAndGetEvent() ; | |
1352 | ||
a1eeb25d | 1353 | if (!fEvent) |
1354 | { | |
6544055e | 1355 | if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry()); |
1356 | return ; | |
1357 | } | |
1358 | ||
3769e0cb | 1359 | //Init pointers, geometry, clusterizer, ocdb, aodb |
6544055e | 1360 | |
3769e0cb | 1361 | InitGeometry(); // only once, must be done before OADB, geo OADB accessed here |
6544055e | 1362 | |
3769e0cb | 1363 | if(fAccessOCDB) AccessOCDB(); |
1364 | if(fAccessOADB) AccessOADB(); // only once | |
1365 | ||
1366 | InitClusterization(); | |
6544055e | 1367 | |
1368 | // Make clusters | |
1369 | if (fJustUnfold) ClusterUnfolding(); | |
1370 | else ClusterizeCells() ; | |
3769e0cb | 1371 | |
6544055e | 1372 | //Recalculate track-matching for the new clusters |
a1eeb25d | 1373 | if(fDoTrackMatching) |
1374 | { | |
1375 | fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom); | |
1376 | } | |
6544055e | 1377 | //Put the new clusters in the AOD list |
1378 | ||
1379 | Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast(); | |
1380 | for(Int_t i = 0; i < kNumberOfCaloClusters; i++){ | |
1381 | AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i); | |
b65291bf | 1382 | |
a1eeb25d | 1383 | newCluster->SetID(i); |
1384 | ||
6544055e | 1385 | //Add matched track |
a1eeb25d | 1386 | if(fDoTrackMatching) |
1387 | { | |
6544055e | 1388 | Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i); |
a1eeb25d | 1389 | if(trackIndex >= 0) |
1390 | { | |
6544055e | 1391 | newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex)); |
1392 | if(DebugLevel() > 1) | |
1393 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i); | |
1394 | } | |
adad4ea9 | 1395 | |
1396 | Float_t dR = 999., dZ = 999.; | |
1397 | fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ); | |
1398 | newCluster->SetTrackDistance(dR,dZ); | |
1399 | ||
6544055e | 1400 | } |
a1eeb25d | 1401 | else |
1402 | {// Assign previously assigned matched track in reco, very very rough | |
adad4ea9 | 1403 | Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster |
1404 | newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]); | |
1405 | } | |
1406 | ||
73a4ed8f | 1407 | //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() ); |
1408 | //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] ); | |
1409 | //printf("\n"); | |
1410 | ||
b96982b6 | 1411 | // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels. |
1412 | fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster); | |
1413 | ||
b65291bf | 1414 | new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster); |
1415 | ||
6544055e | 1416 | if(DebugLevel() > 1 ) |
1417 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E()); | |
1418 | ||
1419 | } // cluster loop | |
1420 | ||
1421 | fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots | |
1422 | ||
1423 | // Clean up | |
1424 | fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why? | |
b65291bf | 1425 | |
6544055e | 1426 | } |
1427 |