]>
Commit | Line | Data |
---|---|---|
eef481c6 | 1 | /************************************************************************** |
6fdd30c4 | 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 | // | |
cd231d42 | 22 | // $Id$ |
6fdd30c4 | 23 | //_________________________________________________________________________ |
24 | ||
25 | // --- Root --- | |
04fcfa08 | 26 | #include <TString.h> |
27 | #include <TRefArray.h> | |
28 | #include <TClonesArray.h> | |
29 | #include <TTree.h> | |
30 | #include <TGeoManager.h> | |
31 | #include <TROOT.h> | |
32 | #include <TInterpreter.h> | |
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" |
6edb4448 | 49 | #include "AliAODMCParticle.h" |
6fdd30c4 | 50 | |
51 | // --- EMCAL | |
6fdd30c4 | 52 | #include "AliEMCALAfterBurnerUF.h" |
53 | #include "AliEMCALGeometry.h" | |
54 | #include "AliEMCALClusterizerNxN.h" | |
55 | #include "AliEMCALClusterizerv1.h" | |
cfcbe5d2 | 56 | #include "AliEMCALClusterizerv2.h" |
6fdd30c4 | 57 | #include "AliEMCALRecPoint.h" |
58 | #include "AliEMCALDigit.h" | |
59 | #include "AliCaloCalibPedestal.h" | |
60 | #include "AliEMCALCalibData.h" | |
61 | ||
62 | #include "AliAnalysisTaskEMCALClusterize.h" | |
63 | ||
64 | ClassImp(AliAnalysisTaskEMCALClusterize) | |
65 | ||
6544055e | 66 | //______________________________________________________________________________ |
6fdd30c4 | 67 | AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name) |
6544055e | 68 | : AliAnalysisTaskSE(name) |
69 | , fEvent(0) | |
3769e0cb | 70 | , fGeom(0), fGeomName("") |
6544055e | 71 | , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE) |
72 | , fCalibData(0), fPedestalData(0) | |
3769e0cb | 73 | , fOCDBpath(""), fAccessOCDB(kFALSE) |
6544055e | 74 | , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0) |
75 | , fRecParam(0), fClusterizer(0) | |
76 | , fUnfolder(0), fJustUnfold(kFALSE) | |
6a797202 | 77 | , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0) |
78 | , fFillAODFile(kFALSE), fFillAODHeader(0) | |
6544055e | 79 | , fFillAODCaloCells(0), fRun(-1) |
80 | , fRecoUtils(0), fConfigName("") | |
6edb4448 | 81 | , fOrgClusterCellId() |
1e0ff9d3 | 82 | , fCellLabels(), fCellSecondLabels(), fCellTime() |
adad4ea9 | 83 | , fCellMatchdEta(), fCellMatchdPhi() |
1e0ff9d3 | 84 | , fRecalibrateWithClusterTime(0) |
3769e0cb | 85 | , fMaxEvent(0), fDoTrackMatching(kFALSE) |
86 | , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0) | |
3e436b3a | 87 | , fRejectBelowThreshold(kFALSE) |
cb68426a | 88 | , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE) |
29d1e4a1 | 89 | , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("") |
cb68426a | 90 | , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("") |
5aa88029 | 91 | , fCentralityClass(""), fSelectEMCALEvent(0) |
92 | , fEMCALEnergyCut(0.), fEMCALNcellsCut (0) | |
6edb4448 | 93 | , fSetCellMCLabelFromCluster(0) |
94 | , fRemapMCLabelForAODs(0) | |
1c197eba | 95 | , fInputFromFilter(0) |
6a8be5c3 | 96 | { |
b43bdd14 | 97 | // Constructor |
98 | ||
e3990982 | 99 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; |
6edb4448 | 100 | |
101 | ResetArrays(); | |
a39f5b70 | 102 | |
cb68426a | 103 | fCentralityBin[0] = fCentralityBin[1]=-1; |
a39f5b70 | 104 | |
6fdd30c4 | 105 | } |
106 | ||
6544055e | 107 | //______________________________________________________________ |
6fdd30c4 | 108 | AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize() |
6544055e | 109 | : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize") |
110 | , fEvent(0) | |
3769e0cb | 111 | , fGeom(0), fGeomName("") |
6544055e | 112 | , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE) |
113 | , fCalibData(0), fPedestalData(0) | |
3769e0cb | 114 | , fOCDBpath(""), fAccessOCDB(kFALSE) |
6544055e | 115 | , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0) |
116 | , fRecParam(0), fClusterizer(0) | |
117 | , fUnfolder(0), fJustUnfold(kFALSE) | |
6a797202 | 118 | , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0) |
119 | , fFillAODFile(kFALSE), fFillAODHeader(0) | |
6544055e | 120 | , fFillAODCaloCells(0), fRun(-1) |
121 | , fRecoUtils(0), fConfigName("") | |
6edb4448 | 122 | , fOrgClusterCellId() |
1e0ff9d3 | 123 | , fCellLabels(), fCellSecondLabels(), fCellTime() |
adad4ea9 | 124 | , fCellMatchdEta(), fCellMatchdPhi() |
1e0ff9d3 | 125 | , fRecalibrateWithClusterTime(0) |
3769e0cb | 126 | , fMaxEvent(0), fDoTrackMatching(kFALSE) |
127 | , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0) | |
3e436b3a | 128 | , fRejectBelowThreshold(kFALSE) |
cb68426a | 129 | , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE) |
29d1e4a1 | 130 | , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("") |
cb68426a | 131 | , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("") |
6c79bfc8 | 132 | , fCentralityClass(""), fSelectEMCALEvent(0) |
133 | , fEMCALEnergyCut(0.), fEMCALNcellsCut (0) | |
6edb4448 | 134 | , fSetCellMCLabelFromCluster(0) |
135 | , fRemapMCLabelForAODs(0) | |
1c197eba | 136 | , fInputFromFilter(0) |
6fdd30c4 | 137 | { |
138 | // Constructor | |
b43bdd14 | 139 | |
e3990982 | 140 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0; |
6edb4448 | 141 | |
142 | ResetArrays(); | |
143 | ||
3af5a372 | 144 | fCentralityBin[0] = fCentralityBin[1]=-1; |
145 | ||
6fdd30c4 | 146 | } |
147 | ||
5994e71f | 148 | |
6544055e | 149 | //_______________________________________________________________ |
6fdd30c4 | 150 | AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize() |
151 | { | |
152 | //dtor | |
153 | ||
b43bdd14 | 154 | if (fDigitsArr) |
155 | { | |
6fdd30c4 | 156 | fDigitsArr->Clear("C"); |
542e1379 | 157 | delete fDigitsArr; |
6fdd30c4 | 158 | } |
159 | ||
b43bdd14 | 160 | if (fClusterArr) |
161 | { | |
6fdd30c4 | 162 | fClusterArr->Delete(); |
542e1379 | 163 | delete fClusterArr; |
6fdd30c4 | 164 | } |
165 | ||
b43bdd14 | 166 | if (fCaloClusterArr) |
167 | { | |
6fdd30c4 | 168 | fCaloClusterArr->Delete(); |
1e0ff9d3 | 169 | delete fCaloClusterArr; |
6fdd30c4 | 170 | } |
6a8be5c3 | 171 | |
385b7abf | 172 | if(fClusterizer) delete fClusterizer; |
173 | if(fUnfolder) delete fUnfolder; | |
174 | if(fRecoUtils) delete fRecoUtils; | |
6a8be5c3 | 175 | |
6fdd30c4 | 176 | } |
177 | ||
5aa88029 | 178 | //_______________________________________________________ |
179 | Bool_t AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() | |
180 | { | |
181 | // Accept event given there is a EMCAL cluster with enough energy, and not noisy, exotic | |
182 | ||
183 | if(!fSelectEMCALEvent) return kTRUE; // accept | |
184 | ||
185 | if(fEMCALEnergyCut <= 0) return kTRUE; // accept | |
186 | ||
f58394f3 | 187 | Int_t nCluster = fEvent -> GetNumberOfCaloClusters(); |
188 | AliVCaloCells * caloCell = fEvent -> GetEMCALCells(); | |
189 | Int_t bc = fEvent -> GetBunchCrossNumber(); | |
190 | ||
5aa88029 | 191 | for(Int_t icalo = 0; icalo < nCluster; icalo++) |
192 | { | |
f58394f3 | 193 | AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo)); |
5aa88029 | 194 | |
195 | if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) && | |
196 | fRecoUtils->IsGoodCluster(clus,fGeom,caloCell,bc)) | |
197 | { | |
198 | ||
199 | if (fDebug > 0) | |
200 | printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Accept : E %2.2f > %2.2f, nCells %d > %d \n", | |
f58394f3 | 201 | clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut); |
5aa88029 | 202 | |
203 | return kTRUE; | |
204 | } | |
205 | ||
206 | }// loop | |
207 | ||
208 | if (fDebug > 0) | |
209 | printf("AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL() - Reject \n"); | |
210 | ||
211 | return kFALSE; | |
212 | ||
f58394f3 | 213 | } |
5aa88029 | 214 | |
3769e0cb | 215 | //_______________________________________________ |
216 | void AliAnalysisTaskEMCALClusterize::AccessOADB() | |
217 | { | |
218 | // Set the AODB calibration, bad channels etc. parameters at least once | |
219 | // alignment matrices from OADB done in SetGeometryMatrices | |
220 | ||
221 | //Set it only once | |
222 | if(fOADBSet) return ; | |
223 | ||
224 | Int_t runnumber = InputEvent()->GetRunNumber() ; | |
225 | TString pass = GetPass(); | |
226 | ||
227 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data()); | |
228 | ||
229 | Int_t nSM = fGeom->GetNumberOfSuperModules(); | |
230 | ||
231 | // Bad map | |
232 | if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()) | |
233 | { | |
234 | AliOADBContainer *contBC=new AliOADBContainer(""); | |
235 | contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels"); | |
236 | ||
237 | TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber); | |
238 | ||
239 | if(arrayBC) | |
240 | { | |
3769e0cb | 241 | fRecoUtils->SwitchOnDistToBadChannelRecalculation(); |
242 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n"); | |
243 | ||
244 | for (Int_t i=0; i<nSM; ++i) | |
245 | { | |
246 | TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i); | |
247 | ||
248 | if (hbm) | |
249 | delete hbm; | |
250 | ||
17e1c80e | 251 | hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i)); |
3769e0cb | 252 | |
253 | if (!hbm) | |
254 | { | |
255 | AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i)); | |
256 | continue; | |
257 | } | |
258 | ||
259 | hbm->SetDirectory(0); | |
260 | fRecoUtils->SetEMCALChannelStatusMap(i,hbm); | |
261 | ||
262 | } // loop | |
17e1c80e | 263 | } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array |
3769e0cb | 264 | } // Remove bad |
265 | ||
266 | // Energy Recalibration | |
267 | if(fRecoUtils->IsRecalibrationOn()) | |
268 | { | |
269 | AliOADBContainer *contRF=new AliOADBContainer(""); | |
270 | ||
271 | contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib"); | |
272 | ||
273 | TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber); | |
274 | ||
275 | if(recal) | |
276 | { | |
277 | TObjArray *recalpass=(TObjArray*)recal->FindObject(pass); | |
278 | ||
279 | if(recalpass) | |
280 | { | |
281 | TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib"); | |
282 | ||
283 | if(recalib) | |
284 | { | |
285 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n"); | |
286 | for (Int_t i=0; i<nSM; ++i) | |
287 | { | |
288 | TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i); | |
289 | ||
290 | if (h) | |
291 | delete h; | |
292 | ||
293 | h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i)); | |
294 | ||
295 | if (!h) | |
296 | { | |
297 | AliError(Form("Could not load EMCALRecalFactors_SM%d",i)); | |
298 | continue; | |
299 | } | |
300 | ||
301 | h->SetDirectory(0); | |
302 | ||
303 | fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h); | |
304 | } // SM loop | |
17e1c80e | 305 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array \n"); // array ok |
306 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok | |
307 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n"); // run number array ok | |
7bf608c9 | 308 | |
3769e0cb | 309 | } // Recalibration on |
310 | ||
7bf608c9 | 311 | // Energy Recalibration, apply on top of previous calibration factors |
312 | if(fRecoUtils->IsRunDepRecalibrationOn()) | |
313 | { | |
314 | AliOADBContainer *contRFTD=new AliOADBContainer(""); | |
315 | ||
316 | contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections"); | |
317 | ||
318 | TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber); | |
319 | ||
74d2310f | 320 | //If it did not exist for this run, get closes one |
321 | if (!htd) | |
322 | { | |
323 | AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber)); | |
324 | // let's get the closest runnumber instead then.. | |
325 | Int_t lower = 0; | |
326 | Int_t ic = 0; | |
327 | Int_t maxEntry = contRFTD->GetNumberOfEntries(); | |
328 | ||
329 | while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) { | |
330 | lower = ic; | |
331 | ic++; | |
332 | } | |
333 | ||
334 | Int_t closest = lower; | |
335 | if ( (ic<maxEntry) && | |
336 | (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) { | |
337 | closest = ic; | |
338 | } | |
339 | ||
340 | AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest))); | |
341 | htd = (TH1S*) contRFTD->GetObjectByIndex(closest); | |
342 | } | |
343 | ||
7bf608c9 | 344 | if(htd) |
345 | { | |
346 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n"); | |
347 | ||
348 | for (Int_t ism=0; ism<nSM; ++ism) | |
349 | { | |
350 | for (Int_t icol=0; icol<48; ++icol) | |
351 | { | |
352 | for (Int_t irow=0; irow<24; ++irow) | |
353 | { | |
354 | Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow); | |
355 | ||
356 | Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor | |
357 | factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T | |
358 | //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID, | |
359 | // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor); | |
360 | fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor); | |
361 | } // columns | |
362 | } // rows | |
363 | } // SM loop | |
364 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n"); | |
365 | } // Run by Run T calibration | |
3769e0cb | 366 | |
367 | // Time Recalibration | |
368 | if(fRecoUtils->IsTimeRecalibrationOn()) | |
369 | { | |
370 | AliOADBContainer *contTRF=new AliOADBContainer(""); | |
371 | ||
f46af216 | 372 | contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib"); |
3769e0cb | 373 | |
374 | TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber); | |
375 | ||
376 | if(trecal) | |
377 | { | |
d4c9aa67 | 378 | TString passM = pass; |
379 | if(pass=="spc_calo") passM = "pass1"; | |
380 | TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM); | |
61841488 | 381 | |
3769e0cb | 382 | if(trecalpass) |
383 | { | |
b43bdd14 | 384 | printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Time Recalibrate EMCAL \n"); |
f46af216 | 385 | for (Int_t ibc = 0; ibc < 4; ++ibc) |
3769e0cb | 386 | { |
f46af216 | 387 | TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc); |
388 | ||
389 | if (h) | |
390 | delete h; | |
391 | ||
392 | h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc)); | |
393 | ||
394 | if (!h) | |
3769e0cb | 395 | { |
f46af216 | 396 | AliError(Form("Could not load hAllTimeAvBC%d",ibc)); |
397 | continue; | |
398 | } | |
399 | ||
400 | h->SetDirectory(0); | |
401 | ||
402 | fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h); | |
403 | } // bunch crossing loop | |
17e1c80e | 404 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok |
405 | }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n"); // run number array ok | |
3769e0cb | 406 | |
407 | } // Time recalibration on | |
408 | ||
409 | // Parameters already set once, so do not it again | |
410 | fOADBSet = kTRUE; | |
411 | ||
412 | } | |
413 | ||
6544055e | 414 | //_________________________________________________ |
415 | Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB() | |
416 | { | |
417 | //Access to OCDB stuff | |
1f77507b | 418 | |
6544055e | 419 | fEvent = InputEvent(); |
420 | if (!fEvent) | |
421 | { | |
422 | Warning("AccessOCDB","Event not available!!!"); | |
423 | return kFALSE; | |
1f77507b | 424 | } |
1f77507b | 425 | |
6544055e | 426 | if (fEvent->GetRunNumber()==fRun) |
427 | return kTRUE; | |
428 | fRun = fEvent->GetRunNumber(); | |
1f77507b | 429 | |
6544055e | 430 | if(DebugLevel() > 1 ) |
596697e6 | 431 | printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Begin"); |
b43bdd14 | 432 | |
6544055e | 433 | AliCDBManager *cdb = AliCDBManager::Instance(); |
1f77507b | 434 | |
1f77507b | 435 | |
b43bdd14 | 436 | if (fOCDBpath.Length()) |
437 | { | |
6544055e | 438 | cdb->SetDefaultStorage(fOCDBpath.Data()); |
439 | printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data()); | |
440 | } | |
441 | ||
442 | cdb->SetRun(fEvent->GetRunNumber()); | |
1f77507b | 443 | |
444 | // | |
6544055e | 445 | // EMCAL from RAW OCDB |
446 | if (fOCDBpath.Contains("alien:")) | |
447 | { | |
448 | cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB"); | |
449 | cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB"); | |
450 | } | |
1f77507b | 451 | |
6544055e | 452 | TString path = cdb->GetDefaultStorage()->GetBaseFolder(); |
1f77507b | 453 | |
6544055e | 454 | // init parameters: |
1f77507b | 455 | |
6544055e | 456 | //Get calibration parameters |
457 | if(!fCalibData) | |
458 | { | |
459 | AliCDBEntry *entry = (AliCDBEntry*) | |
460 | AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); | |
461 | ||
462 | if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); | |
1f77507b | 463 | } |
6544055e | 464 | |
465 | if(!fCalibData) | |
466 | AliFatal("Calibration parameters not found in CDB!"); | |
467 | ||
468 | //Get calibration parameters | |
469 | if(!fPedestalData) | |
470 | { | |
471 | AliCDBEntry *entry = (AliCDBEntry*) | |
472 | AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"); | |
473 | ||
474 | if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject(); | |
1f77507b | 475 | } |
476 | ||
6544055e | 477 | if(!fPedestalData) |
478 | AliFatal("Dead map not found in CDB!"); | |
1f77507b | 479 | |
6544055e | 480 | return kTRUE; |
6fdd30c4 | 481 | } |
482 | ||
9ac44255 | 483 | //_____________________________________________________ |
6544055e | 484 | void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent() |
6fdd30c4 | 485 | { |
6544055e | 486 | // Get the input event, it can depend in embedded events what you want to get |
487 | // Also check if the quality of the event is good if not reject it | |
488 | ||
489 | fEvent = 0x0; | |
f58394f3 | 490 | |
220447cd | 491 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); |
492 | Int_t eventN = Entry(); | |
493 | if(aodIH) eventN = aodIH->GetReadEntry(); | |
494 | ||
6544055e | 495 | if (eventN > fMaxEvent) |
496 | return ; | |
f5036bcb | 497 | |
c50e6949 | 498 | //printf("Clusterizer --- Event %d-- \n",eventN); |
b43bdd14 | 499 | |
6060ed91 | 500 | //Check if input event are embedded events |
501 | //If so, take output event | |
3769e0cb | 502 | if (aodIH && aodIH->GetMergeEvents()) |
503 | { | |
6544055e | 504 | fEvent = AODEvent(); |
1f77507b | 505 | |
506 | if(!aodIH->GetMergeEMCALCells()) | |
507 | AliFatal("Events merged but not EMCAL cells, check analysis settings!"); | |
508 | ||
3769e0cb | 509 | if(DebugLevel() > 1) |
510 | { | |
511 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n"); | |
512 | ||
513 | printf("\t InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(), | |
514 | InputEvent()->GetEMCALCells()->GetNumberOfCells()); | |
515 | ||
516 | printf("\t MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), | |
517 | aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()); | |
518 | ||
519 | for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) | |
520 | { | |
521 | AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl); | |
522 | if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E %f\n",icl,sigCluster->E()); | |
523 | } | |
524 | ||
525 | printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(), | |
526 | AODEvent()->GetEMCALCells()->GetNumberOfCells()); | |
6544055e | 527 | } |
6060ed91 | 528 | } |
1c197eba | 529 | else if(fInputFromFilter) |
530 | { | |
531 | //printf("Get Input From Filtered AOD\n"); | |
532 | fEvent = AODEvent(); | |
533 | } | |
3769e0cb | 534 | else |
535 | { | |
6544055e | 536 | fEvent = InputEvent(); |
1f77507b | 537 | if(fFillAODCaloCells) FillAODCaloCells(); |
538 | if(fFillAODHeader) FillAODHeader(); | |
6060ed91 | 539 | } |
540 | ||
3769e0cb | 541 | if (!fEvent) |
542 | { | |
3fa4fb85 | 543 | Error("UserExec","Event not available"); |
6544055e | 544 | return ; |
6fdd30c4 | 545 | } |
6a8be5c3 | 546 | |
f58394f3 | 547 | //Process events if there is a high energy cluster |
548 | if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; } | |
549 | ||
cd2e4ce6 | 550 | //------------------------------------------------------------------------------------- |
6544055e | 551 | // Reject events if LED was firing, use only for LHC11a data |
552 | // Reject event if triggered by exotic cell and remove exotic cells if not triggered | |
cd2e4ce6 | 553 | //------------------------------------------------------------------------------------- |
6544055e | 554 | |
3769e0cb | 555 | if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; } |
6544055e | 556 | |
3769e0cb | 557 | if( IsExoticEvent() ) { fEvent = 0x0 ; return ; } |
cd2e4ce6 | 558 | |
6a797202 | 559 | //------------------------------------------------------------------------------------- |
560 | // Set the cluster array in the event (output or input) | |
561 | //------------------------------------------------------------------------------------- | |
562 | ||
563 | if ( fFillAODFile ) | |
564 | { | |
565 | //Magic line to write events to AOD filem put after event rejection | |
566 | AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE); | |
567 | } | |
568 | else if( !fOutputAODBranchSet ) | |
569 | { | |
570 | // Create array and put it in the input event, if output AOD not selected, only once | |
571 | InputEvent()->AddObject(fOutputAODBranch); | |
572 | fOutputAODBranchSet = kTRUE; | |
573 | printf("AliAnalysisTaskEMCALClusterize::UserExec() - Add AOD branch <%s> to input event\n",fOutputAODBranchName.Data()); | |
574 | } | |
6544055e | 575 | |
576 | } | |
577 | ||
578 | //____________________________________________________ | |
579 | void AliAnalysisTaskEMCALClusterize::ClusterizeCells() | |
580 | { | |
581 | // Recluster calocells, transform them into digits, | |
582 | // feed the clusterizer with them and get new list of clusters | |
583 | ||
584 | //In case of MC, first loop on the clusters and fill MC label to array | |
585 | Int_t nClusters = fEvent->GetNumberOfCaloClusters(); | |
586 | Int_t nClustersOrg = 0; | |
b43bdd14 | 587 | |
6544055e | 588 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); |
589 | if(aodIH && aodIH->GetEventToMerge()) //Embedding | |
590 | nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal | |
6edb4448 | 591 | |
592 | ResetArrays(); | |
593 | ||
6544055e | 594 | for (Int_t i = 0; i < nClusters; i++) |
595 | { | |
596 | AliVCluster *clus = 0; | |
597 | if(aodIH && aodIH->GetEventToMerge()) //Embedding | |
598 | clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal | |
599 | else | |
600 | clus = fEvent->GetCaloCluster(i); | |
5994e71f | 601 | |
6544055e | 602 | if(!clus) return; |
6a8be5c3 | 603 | |
f46af216 | 604 | if(clus->IsEMCAL()) |
6edb4448 | 605 | { |
6544055e | 606 | Int_t label = clus->GetLabel(); |
607 | Int_t label2 = -1 ; | |
6edb4448 | 608 | //printf("Org cluster E %f, Time %e, Index = %d, ID %d, MC label %d\n", clus->E(), clus->GetTOF(),i, clus->GetID(),label ); |
e59a6a2e | 609 | //printf("Original list of labels from old cluster : \n"); |
6edb4448 | 610 | //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc)); |
e59a6a2e | 611 | |
6544055e | 612 | if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ; |
613 | UShort_t * index = clus->GetCellsAbsId() ; | |
f46af216 | 614 | for(Int_t icell=0; icell < clus->GetNCells(); icell++ ) |
615 | { | |
6edb4448 | 616 | //printf("\t cell %d, MC label %d\n",index[icell],fEvent->GetEMCALCells()->GetCellMCLabel(index[icell])); |
617 | fOrgClusterCellId[index[icell]] = i; | |
6544055e | 618 | fCellLabels[index[icell]] = label; |
619 | fCellSecondLabels[index[icell]] = label2; | |
6edb4448 | 620 | fCellTime[index[icell]] = clus->GetTOF(); |
73a4ed8f | 621 | fCellMatchdEta[index[icell]] = clus->GetTrackDz(); |
622 | fCellMatchdPhi[index[icell]] = clus->GetTrackDx(); | |
6fdd30c4 | 623 | } |
6544055e | 624 | nClustersOrg++; |
6fdd30c4 | 625 | } |
b43bdd14 | 626 | // printf("\n"); |
6544055e | 627 | } |
628 | ||
629 | // Transform CaloCells into Digits | |
b43bdd14 | 630 | |
6544055e | 631 | Int_t idigit = 0; |
632 | Int_t id = -1; | |
633 | Float_t amp = -1; | |
634 | Double_t time = -1; | |
635 | ||
6edb4448 | 636 | AliVCaloCells *cells = fEvent->GetEMCALCells(); |
b43bdd14 | 637 | |
a7e5a381 | 638 | Int_t bc = InputEvent()->GetBunchCrossNumber(); |
11e7733f | 639 | |
6544055e | 640 | for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) |
641 | { | |
6544055e | 642 | // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell |
643 | id = cells->GetCellNumber(icell); | |
a7e5a381 | 644 | Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells); |
645 | ||
1e0ff9d3 | 646 | // Do not include cells with too low energy, nor exotic cell |
d95102d2 | 647 | if( amp < fRecParam->GetMinECut() || |
648 | time > fRecParam->GetTimeMax() || | |
649 | time < fRecParam->GetTimeMin() ) accept = kFALSE; | |
1e0ff9d3 | 650 | |
651 | // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs. | |
652 | if (fRecalibrateWithClusterTime) | |
653 | { | |
654 | time = fCellTime[id]; | |
655 | //printf("cell %d time org %f - ",id, time*1.e9); | |
656 | fRecoUtils->RecalibrateCellTime(id,bc,time); | |
657 | //printf("recal %f\n",time*1.e9); | |
658 | } | |
a7e5a381 | 659 | |
d95102d2 | 660 | //Exotic? |
661 | if (accept && fRecoUtils->IsExoticCell(id,cells,bc)) | |
662 | accept = kFALSE; | |
b43bdd14 | 663 | |
f46af216 | 664 | if( !accept ) |
665 | { | |
f46af216 | 666 | if( DebugLevel() > 2 ) |
667 | printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n", | |
b43bdd14 | 668 | id,icell, cells->GetNumberOfCells(), amp, time*1.e9); |
6544055e | 669 | continue; |
670 | } | |
b43bdd14 | 671 | |
bdd50327 | 672 | Int_t mcLabel = cells->GetMCLabel(icell); |
e59a6a2e | 673 | //printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - cell %d, mc label %d\n",id,mcLabel); |
6edb4448 | 674 | |
bdd50327 | 675 | //if(fCellLabels[id]!=mcLabel)printf("mcLabel %d - %d\n",mcLabel,fCellLabels[id]); |
6edb4448 | 676 | if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions |
677 | else if( fSetCellMCLabelFromCluster == 0 && fRemapMCLabelForAODs) RemapMCLabelForAODs(mcLabel); | |
e59a6a2e | 678 | else mcLabel = -1; // found later |
679 | ||
680 | //printf("\t new label %d\n",mcLabel); | |
bdd50327 | 681 | |
682 | // Create the digit, put a fake primary deposited energy to trick the clusterizer | |
683 | // when checking the most likely primary | |
11e7733f | 684 | |
685 | Float_t efrac = cells->GetEFraction(icell); | |
686 | ||
687 | //When checking the MC of digits, give weight to cells with embedded signal | |
688 | if (mcLabel > 0 && efrac < 1.e-6) efrac = 1; | |
689 | ||
690 | //printf("******* Cell %d, id %d, e %f, fraction %f, MC label %d, used MC label %d\n",icell,id,amp,cells->GetEFraction(icell),cells->GetMCLabel(icell),mcLabel); | |
691 | ||
692 | new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, amp*efrac); | |
e59a6a2e | 693 | // Last parameter should be MC deposited energy, since it is not available, add just the cell amplitude so that |
694 | // we give more weight to the MC label of the cell with highest energy in the cluster | |
6edb4448 | 695 | |
6544055e | 696 | idigit++; |
697 | } | |
6a8be5c3 | 698 | |
04fcfa08 | 699 | fDigitsArr->Sort(); |
b43bdd14 | 700 | |
6544055e | 701 | //------------------------------------------------------------------------------------- |
702 | //Do the clusterization | |
703 | //------------------------------------------------------------------------------------- | |
b43bdd14 | 704 | |
6544055e | 705 | fClusterizer->Digits2Clusters(""); |
706 | ||
707 | //------------------------------------------------------------------------------------- | |
708 | //Transform the recpoints into AliVClusters | |
709 | //------------------------------------------------------------------------------------- | |
710 | ||
b43bdd14 | 711 | RecPoints2Clusters(); |
6544055e | 712 | |
3769e0cb | 713 | if(!fCaloClusterArr) |
714 | { | |
6544055e | 715 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast()); |
716 | return; | |
717 | } | |
718 | ||
3769e0cb | 719 | if( DebugLevel() > 0 ) |
720 | { | |
6544055e | 721 | printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast()); |
b43bdd14 | 722 | |
723 | if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()) | |
724 | { | |
6544055e | 725 | printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n", |
726 | fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(), | |
727 | fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries()); | |
728 | } | |
729 | } | |
6544055e | 730 | } |
a39f5b70 | 731 | |
73a4ed8f | 732 | //_____________________________________________________ |
6544055e | 733 | void AliAnalysisTaskEMCALClusterize::ClusterUnfolding() |
734 | { | |
735 | // Take the event clusters and unfold them | |
736 | ||
737 | AliVCaloCells *cells = fEvent->GetEMCALCells(); | |
738 | Double_t cellAmplitude = 0; | |
739 | Double_t cellTime = 0; | |
740 | Short_t cellNumber = 0; | |
48d9e9e7 | 741 | Int_t cellMCLabel = 0; |
77e93dc2 | 742 | Double_t cellEFrac = 0; |
6544055e | 743 | Int_t nClustersOrg = 0; |
744 | ||
745 | // Fill the array with the EMCAL clusters, copy them | |
746 | for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++) | |
747 | { | |
748 | AliVCluster *clus = fEvent->GetCaloCluster(i); | |
b43bdd14 | 749 | if(clus->IsEMCAL()) |
750 | { | |
6544055e | 751 | //recalibrate/remove bad channels/etc if requested |
b43bdd14 | 752 | if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())) |
753 | { | |
5994e71f | 754 | continue; |
6544055e | 755 | } |
3bfc4732 | 756 | |
b43bdd14 | 757 | if(fRecoUtils->IsRecalibrationOn()) |
758 | { | |
6544055e | 759 | //Calibrate cluster |
760 | fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells); | |
761 | ||
762 | //CalibrateCells | |
763 | for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) | |
764 | { | |
77e93dc2 | 765 | if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE) |
6544055e | 766 | break; |
767 | ||
768 | Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; | |
769 | fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta); | |
770 | fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); | |
771 | ||
772 | //Do not include bad channels found in analysis? | |
773 | if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() && | |
6edb4448 | 774 | fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)) |
6544055e | 775 | continue; |
6544055e | 776 | |
777 | cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime); | |
778 | ||
779 | }// cells loop | |
780 | }// recalibrate | |
5994e71f | 781 | |
6544055e | 782 | //Cast to ESD or AOD, needed to create the cluster array |
783 | AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus); | |
784 | AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus); | |
b43bdd14 | 785 | |
786 | if (esdCluster) | |
787 | { | |
6544055e | 788 | fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) ); |
789 | }//ESD | |
b43bdd14 | 790 | else if(aodCluster) |
791 | { | |
6544055e | 792 | fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) ); |
793 | }//AOD | |
794 | else | |
795 | Warning("UserExec()"," - Wrong CaloCluster type?"); | |
b43bdd14 | 796 | |
6544055e | 797 | nClustersOrg++; |
6a8be5c3 | 798 | } |
6fdd30c4 | 799 | } |
800 | ||
6544055e | 801 | //Do the unfolding |
802 | fUnfolder->UnfoldClusters(fCaloClusterArr, cells); | |
6a8be5c3 | 803 | |
6544055e | 804 | //CLEAN-UP |
805 | fUnfolder->Clear(); | |
385b7abf | 806 | |
6544055e | 807 | } |
808 | ||
809 | //_____________________________________________________ | |
810 | void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() | |
811 | { | |
812 | // Put calo cells in standard branch | |
813 | AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells()); | |
814 | Int_t nEMcell = eventEMcells.GetNumberOfCells() ; | |
6fdd30c4 | 815 | |
6544055e | 816 | AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells()); |
817 | aodEMcells.CreateContainer(nEMcell); | |
818 | aodEMcells.SetType(AliVCaloCells::kEMCALCell); | |
819 | Double_t calibFactor = 1.; | |
b43bdd14 | 820 | for (Int_t iCell = 0; iCell < nEMcell; iCell++) |
821 | { | |
6544055e | 822 | Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1; |
823 | fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta); | |
824 | fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta); | |
385b7abf | 825 | |
b43bdd14 | 826 | if(fRecoUtils->IsRecalibrationOn()) |
827 | { | |
6544055e | 828 | calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi); |
385b7abf | 829 | } |
0821ba9f | 830 | |
b43bdd14 | 831 | if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)) |
832 | { //Channel is not declared as bad | |
77e93dc2 | 833 | aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor, |
834 | eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell)); | |
6544055e | 835 | } |
b43bdd14 | 836 | else |
837 | { | |
77e93dc2 | 838 | aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0); |
5994e71f | 839 | } |
6fdd30c4 | 840 | } |
6544055e | 841 | aodEMcells.Sort(); |
6fdd30c4 | 842 | |
6544055e | 843 | } |
6fdd30c4 | 844 | |
6544055e | 845 | //__________________________________________________ |
846 | void AliAnalysisTaskEMCALClusterize::FillAODHeader() | |
6fdd30c4 | 847 | { |
6544055e | 848 | //Put event header information in standard AOD branch |
e9dd2d80 | 849 | |
6544055e | 850 | AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent); |
851 | AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent); | |
e9dd2d80 | 852 | |
6544055e | 853 | Double_t pos[3] ; |
854 | Double_t covVtx[6]; | |
855 | for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.; | |
856 | ||
857 | AliAODHeader* header = AODEvent()->GetHeader(); | |
858 | header->SetRunNumber(fEvent->GetRunNumber()); | |
859 | ||
b43bdd14 | 860 | if(esdevent) |
861 | { | |
6544055e | 862 | TTree* tree = fInputHandler->GetTree(); |
b43bdd14 | 863 | if (tree) |
864 | { | |
6544055e | 865 | TFile* file = tree->GetCurrentFile(); |
866 | if (file) header->SetESDFileName(file->GetName()); | |
867 | } | |
868 | } | |
869 | else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName()); | |
e9dd2d80 | 870 | |
6544055e | 871 | header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber()); |
872 | header->SetOrbitNumber(fEvent->GetOrbitNumber()); | |
873 | header->SetPeriodNumber(fEvent->GetPeriodNumber()); | |
874 | header->SetEventType(fEvent->GetEventType()); | |
385b7abf | 875 | |
6544055e | 876 | //Centrality |
b43bdd14 | 877 | if(fEvent->GetCentrality()) |
6544055e | 878 | header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality()))); |
b43bdd14 | 879 | else |
6544055e | 880 | header->SetCentrality(0); |
385b7abf | 881 | |
6544055e | 882 | //Trigger |
883 | header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection | |
997b261e | 884 | |
885 | header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses()); | |
b43bdd14 | 886 | |
6544055e | 887 | header->SetTriggerMask(fEvent->GetTriggerMask()); |
888 | header->SetTriggerCluster(fEvent->GetTriggerCluster()); | |
b43bdd14 | 889 | |
890 | if (esdevent) | |
891 | { | |
6544055e | 892 | header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs()); |
893 | header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs()); | |
894 | header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs()); | |
895 | } | |
b43bdd14 | 896 | else if (aodevent) |
897 | { | |
6544055e | 898 | header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs()); |
899 | header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs()); | |
900 | header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs()); | |
901 | } | |
e9dd2d80 | 902 | |
6544055e | 903 | header->SetMagneticField(fEvent->GetMagneticField()); |
904 | //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.); | |
e9dd2d80 | 905 | |
6544055e | 906 | header->SetZDCN1Energy(fEvent->GetZDCN1Energy()); |
907 | header->SetZDCP1Energy(fEvent->GetZDCP1Energy()); | |
908 | header->SetZDCN2Energy(fEvent->GetZDCN2Energy()); | |
909 | header->SetZDCP2Energy(fEvent->GetZDCP2Energy()); | |
910 | header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1)); | |
385b7abf | 911 | |
76705ce7 | 912 | Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()}; |
6544055e | 913 | Float_t diamcov[3]; |
914 | fEvent->GetDiamondCovXY(diamcov); | |
915 | header->SetDiamond(diamxy,diamcov); | |
916 | if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ()); | |
917 | else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ()); | |
e9dd2d80 | 918 | |
6544055e | 919 | // |
920 | Int_t nVertices = 1 ;/* = prim. vtx*/; | |
921 | Int_t nCaloClus = fEvent->GetNumberOfCaloClusters(); | |
e9dd2d80 | 922 | |
6544055e | 923 | AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0); |
385b7abf | 924 | |
6544055e | 925 | // Access to the AOD container of vertices |
926 | TClonesArray &vertices = *(AODEvent()->GetVertices()); | |
927 | Int_t jVertices=0; | |
385b7abf | 928 | |
6544055e | 929 | // Add primary vertex. The primary tracks will be defined |
930 | // after the loops on the composite objects (V0, cascades, kinks) | |
931 | fEvent->GetPrimaryVertex()->GetXYZ(pos); | |
932 | Float_t chi = 0; | |
b43bdd14 | 933 | if (esdevent) |
934 | { | |
6544055e | 935 | esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx); |
936 | chi = esdevent->GetPrimaryVertex()->GetChi2toNDF(); | |
937 | } | |
b43bdd14 | 938 | else if (aodevent) |
939 | { | |
6544055e | 940 | aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx); |
941 | chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD? | |
6fdd30c4 | 942 | } |
bd9e8ebd | 943 | |
6544055e | 944 | AliAODVertex * primary = new(vertices[jVertices++]) |
945 | AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary); | |
946 | primary->SetName(fEvent->GetPrimaryVertex()->GetName()); | |
947 | primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle()); | |
e9dd2d80 | 948 | |
6544055e | 949 | } |
950 | ||
afc22931 | 951 | //___________________________________________________________ |
952 | void AliAnalysisTaskEMCALClusterize::FillCaloClusterInEvent() | |
b43bdd14 | 953 | { |
954 | // Get the CaloClusters array, do some final calculations | |
afc22931 | 955 | // and put the clusters in the output or input event |
956 | // as a separate branch | |
b43bdd14 | 957 | |
958 | //First recalculate track-matching for the new clusters | |
959 | if(fDoTrackMatching) | |
960 | { | |
961 | fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom); | |
962 | } | |
963 | //Put the new clusters in the AOD list | |
964 | ||
965 | Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast(); | |
596697e6 | 966 | |
b43bdd14 | 967 | for(Int_t i = 0; i < kNumberOfCaloClusters; i++) |
968 | { | |
969 | AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i); | |
970 | ||
971 | newCluster->SetID(i); | |
972 | ||
afc22931 | 973 | // Correct cluster energy non linearity |
974 | newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster)); | |
975 | ||
b43bdd14 | 976 | //Add matched track |
977 | if(fDoTrackMatching) | |
978 | { | |
979 | Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i); | |
980 | if(trackIndex >= 0) | |
981 | { | |
982 | newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex)); | |
983 | if(DebugLevel() > 1) | |
984 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i); | |
985 | } | |
986 | ||
987 | Float_t dR = 999., dZ = 999.; | |
988 | fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ); | |
989 | newCluster->SetTrackDistance(dR,dZ); | |
990 | ||
991 | } | |
992 | else | |
993 | {// Assign previously assigned matched track in reco, very very rough | |
994 | Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster | |
995 | newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]); | |
996 | } | |
1e0ff9d3 | 997 | |
998 | //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() ); | |
999 | //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] ); | |
1000 | //printf("\n"); | |
1001 | ||
b43bdd14 | 1002 | // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels. |
1003 | fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster); | |
1004 | ||
1005 | new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster); | |
1006 | ||
11e7733f | 1007 | if(DebugLevel() > 1 ) |
1008 | printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f, mc label %d \n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()); | |
b43bdd14 | 1009 | |
1010 | } // cluster loop | |
1011 | ||
1012 | fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots | |
1013 | ||
b43bdd14 | 1014 | |
1015 | } | |
1016 | ||
3769e0cb | 1017 | //_______________________________________________ |
1018 | TString AliAnalysisTaskEMCALClusterize::GetPass() | |
1019 | { | |
1020 | // Get passx from filename. | |
1021 | ||
1022 | if (!AliAnalysisManager::GetAnalysisManager()->GetTree()) | |
1023 | { | |
1024 | AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n"); | |
1025 | return TString(""); | |
1026 | } | |
1027 | ||
1028 | if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()) | |
1029 | { | |
1030 | AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n"); | |
1031 | return TString(""); | |
1032 | } | |
1033 | ||
1034 | TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName()); | |
1035 | if (pass.Contains("ass1")) return TString("pass1"); | |
1036 | else if (pass.Contains("ass2")) return TString("pass2"); | |
1037 | else if (pass.Contains("ass3")) return TString("pass3"); | |
61841488 | 1038 | else if (pass.Contains("ass4")) return TString("pass4"); |
1039 | else if (pass.Contains("ass5")) return TString("pass5"); | |
9783326d | 1040 | else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo"); |
a4976005 | 1041 | else if (pass.Contains("calo") || pass.Contains("high_lumi")) |
1042 | { | |
1043 | printf("AliAnalysisTaskEMCALClusterize::GetPass() - Path contains <calo> or <high-lumi>, set as <pass1>\n"); | |
1044 | return TString("pass1"); | |
1045 | } | |
3769e0cb | 1046 | |
1047 | // No condition fullfilled, give a default value | |
1048 | printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n"); | |
1049 | return TString(""); | |
1050 | ||
1051 | } | |
1052 | ||
6544055e | 1053 | //_________________________________________ |
1054 | void AliAnalysisTaskEMCALClusterize::Init() | |
1055 | { | |
1056 | //Init analysis with configuration macro if available | |
1057 | ||
3769e0cb | 1058 | fOADBSet = kFALSE; |
1059 | if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ; | |
1060 | ||
3769e0cb | 1061 | fBranchNames = "ESD:AliESDHeader.,EMCALCells."; |
1062 | ||
1063 | if(!fRecParam) fRecParam = new AliEMCALRecParam; | |
1064 | if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils(); | |
1065 | ||
1066 | if(fMaxEvent <= 0) fMaxEvent = 1000000000; | |
1067 | if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005; | |
1068 | if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001; | |
3e436b3a | 1069 | fRejectBelowThreshold = kFALSE; |
1070 | ||
cb68426a | 1071 | //Centrality |
1072 | if(fCentralityClass == "") fCentralityClass = "V0M"; | |
cb68426a | 1073 | |
3769e0cb | 1074 | if (fOCDBpath == "") fOCDBpath = "raw://" ; |
1075 | if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ; | |
1076 | ||
1077 | if(gROOT->LoadMacro(fConfigName) >=0) | |
1078 | { | |
6544055e | 1079 | printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data()); |
1080 | AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()"); | |
1081 | fGeomName = clus->fGeomName; | |
1082 | fLoadGeomMatrices = clus->fLoadGeomMatrices; | |
1083 | fOCDBpath = clus->fOCDBpath; | |
1084 | fAccessOCDB = clus->fAccessOCDB; | |
1085 | fRecParam = clus->fRecParam; | |
1086 | fJustUnfold = clus->fJustUnfold; | |
1087 | fFillAODFile = clus->fFillAODFile; | |
1088 | fRecoUtils = clus->fRecoUtils; | |
1089 | fConfigName = clus->fConfigName; | |
1090 | fMaxEvent = clus->fMaxEvent; | |
1091 | fDoTrackMatching = clus->fDoTrackMatching; | |
1092 | fOutputAODBranchName = clus->fOutputAODBranchName; | |
e3990982 | 1093 | for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ; |
cb68426a | 1094 | fCentralityClass = clus->fCentralityClass; |
1095 | fCentralityBin[0] = clus->fCentralityBin[0]; | |
1096 | fCentralityBin[1] = clus->fCentralityBin[1]; | |
6fdd30c4 | 1097 | } |
29d1e4a1 | 1098 | |
6544055e | 1099 | } |
6fdd30c4 | 1100 | |
6544055e | 1101 | //_______________________________________________________ |
6fdd30c4 | 1102 | void AliAnalysisTaskEMCALClusterize::InitClusterization() |
1103 | { | |
1104 | //Select clusterization/unfolding algorithm and set all the needed parameters | |
1105 | ||
04fcfa08 | 1106 | if (fJustUnfold) |
1107 | { | |
bd9e8ebd | 1108 | // init the unfolding afterburner |
1109 | delete fUnfolder; | |
1759f743 | 1110 | fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut()); |
bd9e8ebd | 1111 | return; |
6a8be5c3 | 1112 | } |
1113 | ||
bd9e8ebd | 1114 | //First init the clusterizer |
1115 | delete fClusterizer; | |
1116 | if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1) | |
1117 | fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData); | |
cfcbe5d2 | 1118 | else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2) |
1119 | fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData); | |
b43bdd14 | 1120 | else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) |
1121 | { | |
bd9e8ebd | 1122 | fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData); |
212ac797 | 1123 | fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); |
1124 | fClusterizer->SetNColDiff(fRecParam->GetNColDiff()); | |
04fcfa08 | 1125 | } |
1126 | else | |
1127 | { | |
bd9e8ebd | 1128 | AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag())); |
1129 | } | |
6a8be5c3 | 1130 | |
bd9e8ebd | 1131 | //Now set the parameters |
1132 | fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() ); | |
1133 | fClusterizer->SetECALogWeight ( fRecParam->GetW0() ); | |
1134 | fClusterizer->SetMinECut ( fRecParam->GetMinECut() ); | |
1135 | fClusterizer->SetUnfolding ( fRecParam->GetUnfold() ); | |
1136 | fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() ); | |
1137 | fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() ); | |
1138 | fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() ); | |
1139 | fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() ); | |
1140 | fClusterizer->SetInputCalibrated ( kTRUE ); | |
2af35150 | 1141 | fClusterizer->SetJustClusters ( kTRUE ); |
04fcfa08 | 1142 | |
1143 | // Initialize the cluster rec points and digits arrays and get them. | |
1144 | fClusterizer->SetOutput(0); | |
1145 | fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints()); | |
1146 | fDigitsArr = fClusterizer->GetDigits(); | |
1147 | ||
bd9e8ebd | 1148 | //In case of unfolding after clusterization is requested, set the corresponding parameters |
04fcfa08 | 1149 | if(fRecParam->GetUnfold()) |
1150 | { | |
bd9e8ebd | 1151 | Int_t i=0; |
04fcfa08 | 1152 | for (i = 0; i < 8; i++) |
1153 | { | |
bd9e8ebd | 1154 | fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i)); |
1155 | }//end of loop over parameters | |
04fcfa08 | 1156 | |
1157 | for (i = 0; i < 3; i++) | |
1158 | { | |
bd9e8ebd | 1159 | fClusterizer->SetPar5 (i, fRecParam->GetPar5(i)); |
1160 | fClusterizer->SetPar6 (i, fRecParam->GetPar6(i)); | |
1161 | }//end of loop over parameters | |
3e436b3a | 1162 | fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy |
bd9e8ebd | 1163 | fClusterizer->InitClusterUnfolding(); |
6544055e | 1164 | |
bd9e8ebd | 1165 | }// to unfold |
6fdd30c4 | 1166 | } |
3fa4fb85 | 1167 | |
6544055e | 1168 | //_________________________________________________ |
1169 | void AliAnalysisTaskEMCALClusterize::InitGeometry() | |
1170 | { | |
3769e0cb | 1171 | // Init geometry and set the geometry matrix, for the first event, skip the rest |
6544055e | 1172 | // Also set once the run dependent calibrations |
1173 | ||
3769e0cb | 1174 | if(fGeomMatrixSet) return; |
1175 | ||
1176 | Int_t runnumber = InputEvent()->GetRunNumber() ; | |
1177 | if (!fGeom) | |
1178 | { | |
1179 | if(fGeomName=="") | |
1180 | { | |
1181 | if (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1"; | |
1182 | else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1"; | |
1183 | else fGeomName = "EMCAL_COMPLETE12SMV1"; | |
29d1e4a1 | 1184 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n", |
1185 | fGeomName.Data(),runnumber); | |
3769e0cb | 1186 | } |
6544055e | 1187 | |
3769e0cb | 1188 | fGeom = AliEMCALGeometry::GetInstance(fGeomName); |
6544055e | 1189 | |
29d1e4a1 | 1190 | // Init geometry, I do not like much to do it like this ... |
1191 | if(fImportGeometryFromFile && !gGeoManager) | |
1192 | { | |
1193 | if(fImportGeometryFilePath=="") // If not specified, set location depending on run number | |
1194 | { | |
1195 | // "$ALICE_ROOT/EVE/alice-data/default_geo.root" | |
1196 | if (runnumber < 140000 && | |
1197 | runnumber >= 100000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2010.root"; | |
1198 | if (runnumber >= 140000 && | |
1199 | runnumber < 171000) fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2011.root"; | |
1200 | else fImportGeometryFilePath = "$ALICE_ROOT/OADB/EMCAL/geometry_2012.root"; // 2012-2013 | |
1201 | } | |
1202 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Import %s\n",fImportGeometryFilePath.Data()); | |
1203 | TGeoManager::Import(fImportGeometryFilePath) ; | |
1204 | } | |
1205 | ||
3769e0cb | 1206 | if(fDebug > 0) |
1207 | { | |
1208 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber); | |
1209 | if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices"); | |
1210 | printf("\n"); | |
1211 | } | |
1212 | } // geometry pointer did not exist before | |
1213 | ||
1214 | if(fLoadGeomMatrices) | |
1215 | { | |
1216 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n"); | |
1217 | ||
1218 | // OADB if available | |
1219 | AliOADBContainer emcGeoMat("AliEMCALgeo"); | |
1220 | emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo"); | |
1221 | TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices"); | |
1222 | ||
1223 | ||
1224 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++) | |
1225 | { | |
1226 | ||
1227 | if (!fGeomMatrix[mod]) // Get it from OADB | |
1228 | { | |
1229 | if(fDebug > 1 ) | |
1230 | printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n", | |
1231 | mod,((TGeoHMatrix*) matEMCAL->At(mod))); | |
1232 | //((TGeoHMatrix*) matEMCAL->At(mod))->Print(); | |
1233 | ||
1234 | fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ; | |
1235 | } | |
1236 | ||
1237 | if(fGeomMatrix[mod]) | |
1238 | { | |
1239 | if(DebugLevel() > 1) | |
1240 | fGeomMatrix[mod]->Print(); | |
1241 | ||
1242 | fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ; | |
1243 | } | |
1244 | ||
1245 | fGeomMatrixSet=kTRUE; | |
1246 | ||
1247 | }//SM loop | |
1248 | }//Load matrices | |
1249 | else if(!gGeoManager) | |
1250 | { | |
1251 | printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data"); | |
1252 | //Still not implemented in AOD, just a workaround to be able to work at least with ESDs | |
1253 | if(!strcmp(fEvent->GetName(),"AliAODEvent")) | |
1254 | { | |
1255 | if(DebugLevel() > 1) | |
1256 | Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs."); | |
1257 | }//AOD | |
1258 | else | |
1259 | { | |
1260 | if(DebugLevel() > 1) | |
1261 | printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices."); | |
1262 | ||
1263 | AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ; | |
1264 | ||
1265 | if(!esd) | |
1266 | { | |
1267 | Error("InitGeometry"," - This event does not contain ESDs?"); | |
6a797202 | 1268 | if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE); |
3769e0cb | 1269 | return; |
1270 | } | |
1271 | ||
1272 | for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++) | |
1273 | { | |
1274 | if(DebugLevel() > 1) | |
1275 | esd->GetEMCALMatrix(mod)->Print(); | |
1276 | ||
1277 | if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ; | |
1278 | ||
1279 | } | |
1280 | ||
1281 | fGeomMatrixSet=kTRUE; | |
1282 | ||
1283 | }//ESD | |
1284 | }//Load matrices from Data | |
1285 | ||
6544055e | 1286 | } |
1287 | ||
6544055e | 1288 | //____________________________________________________ |
1289 | Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent() | |
1290 | { | |
1291 | ||
1292 | // Check if event is exotic, get an exotic cell and compare with triggered patch | |
1293 | // If there is a match remove event ... to be completed, filled with something provisional | |
1294 | ||
1295 | if(!fRemoveExoticEvents) return kFALSE; | |
1296 | ||
1297 | // Loop on cells | |
1298 | AliVCaloCells * cells = fEvent->GetEMCALCells(); | |
1299 | Float_t totCellE = 0; | |
a7e5a381 | 1300 | Int_t bc = InputEvent()->GetBunchCrossNumber(); |
b43bdd14 | 1301 | for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) |
1302 | { | |
6544055e | 1303 | |
1304 | Float_t ecell = 0 ; | |
1305 | Double_t tcell = 0 ; | |
1306 | ||
1307 | Int_t absID = cells->GetCellNumber(icell); | |
a7e5a381 | 1308 | Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells); |
1309 | if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell; | |
6544055e | 1310 | } |
1311 | ||
997b261e | 1312 | // TString triggerclasses = event->GetFiredTriggerClasses(); |
6544055e | 1313 | // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4); |
6a797202 | 1314 | // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);; |
6544055e | 1315 | // return; |
1316 | // | |
1317 | ||
1318 | //printf("TotE cell %f\n",totCellE); | |
1319 | if(totCellE < 1) return kTRUE; | |
1320 | ||
1321 | return kFALSE; | |
1322 | ||
1323 | } | |
1324 | ||
3769e0cb | 1325 | //________________________________________________________________ |
1326 | Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run) | |
6544055e | 1327 | { |
1328 | //Check if event is LED | |
1329 | ||
1330 | if(!fRemoveLEDEvents) return kFALSE; | |
3769e0cb | 1331 | |
1332 | //check events of LHC11a period | |
b9d6234e | 1333 | if(run < 146858 || run > 146860) return kFALSE ; |
3769e0cb | 1334 | |
6544055e | 1335 | // Count number of cells with energy larger than 0.1 in SM3, cut on this number |
1336 | Int_t ncellsSM3 = 0; | |
6544055e | 1337 | AliVCaloCells * cells = fEvent->GetEMCALCells(); |
3769e0cb | 1338 | for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++) |
1339 | { | |
6544055e | 1340 | if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++; |
6544055e | 1341 | } |
1342 | ||
997b261e | 1343 | TString triggerclasses = fEvent->GetFiredTriggerClasses(); |
6544055e | 1344 | |
1345 | Int_t ncellcut = 21; | |
1346 | if(triggerclasses.Contains("EMC")) ncellcut = 35; | |
1347 | ||
b9d6234e | 1348 | if( ncellsSM3 >= ncellcut) |
3769e0cb | 1349 | { |
b9d6234e | 1350 | printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d\n",(Int_t)Entry(),ncellsSM3); |
6a797202 | 1351 | if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);; |
6544055e | 1352 | return kTRUE; |
1353 | } | |
1354 | ||
1355 | return kFALSE; | |
1356 | ||
1357 | } | |
1358 | ||
b43bdd14 | 1359 | //_______________________________________________________ |
1360 | void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() | |
6fdd30c4 | 1361 | { |
1362 | // Restore clusters from recPoints | |
b43bdd14 | 1363 | // Cluster energy, global position, cells and their amplitude |
1364 | // fractions are restored | |
1365 | ||
6a8be5c3 | 1366 | Int_t j = 0; |
b43bdd14 | 1367 | for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++) |
6fdd30c4 | 1368 | { |
b43bdd14 | 1369 | AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i); |
6fdd30c4 | 1370 | |
1371 | const Int_t ncells = recPoint->GetMultiplicity(); | |
f6e45fe0 | 1372 | Int_t ncellsTrue = 0; |
6fdd30c4 | 1373 | |
6a8be5c3 | 1374 | if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue; |
1375 | ||
6fdd30c4 | 1376 | // cells and their amplitude fractions |
1377 | UShort_t absIds[ncells]; | |
1378 | Double32_t ratios[ncells]; | |
1379 | ||
f6e45fe0 | 1380 | //For later check embedding |
1381 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); | |
f053064e | 1382 | |
1383 | Float_t clusterE = 0; | |
f46af216 | 1384 | for (Int_t c = 0; c < ncells; c++) |
1385 | { | |
b43bdd14 | 1386 | AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]); |
6fdd30c4 | 1387 | |
f6e45fe0 | 1388 | absIds[ncellsTrue] = digit->GetId(); |
1389 | ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude(); | |
11e7733f | 1390 | |
d390e7cb | 1391 | // In case of unfolding, remove digits with unfolded energy too low |
b43bdd14 | 1392 | if(fSelectCell) |
1393 | { | |
f46af216 | 1394 | if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac) |
1395 | { | |
11e7733f | 1396 | if(DebugLevel() > 1) |
f46af216 | 1397 | { |
6544055e | 1398 | printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n", |
d390e7cb | 1399 | recPoint->GetEnergiesList()[c],digit->GetAmplitude()); |
1400 | } | |
1401 | ||
1402 | continue; | |
1403 | ||
1404 | } // if cuts | |
1405 | }// Select cells | |
6544055e | 1406 | |
d390e7cb | 1407 | //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected |
d390e7cb | 1408 | clusterE +=recPoint->GetEnergiesList()[c]; |
f053064e | 1409 | |
f6e45fe0 | 1410 | // In case of embedding, fill ratio with amount of signal, |
f46af216 | 1411 | if (aodIH && aodIH->GetMergeEvents()) |
1412 | { | |
f6e45fe0 | 1413 | //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells(); |
1414 | AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells(); | |
1415 | AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells(); | |
1416 | ||
1417 | Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]); | |
1418 | //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]); | |
1419 | Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]); | |
1420 | //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude()); | |
1421 | ||
1422 | if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude; | |
1423 | //printf("\t \t ratio %f\n",ratios[ncellsTrue]); | |
1424 | ||
1425 | }//Embedding | |
1426 | ||
11e7733f | 1427 | ncellsTrue++; |
1428 | ||
f6e45fe0 | 1429 | }// cluster cell loop |
6fdd30c4 | 1430 | |
f46af216 | 1431 | if (ncellsTrue < 1) |
1432 | { | |
f053064e | 1433 | if (DebugLevel() > 1) |
1434 | printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n", | |
1435 | recPoint->GetEnergy(), ncells); | |
1436 | continue; | |
1437 | } | |
1438 | ||
1439 | //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue); | |
1440 | ||
f46af216 | 1441 | if(clusterE < fRecParam->GetClusteringThreshold()) |
1442 | { | |
f053064e | 1443 | if (DebugLevel()>1) |
1444 | printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE); | |
6fdd30c4 | 1445 | continue; |
1446 | } | |
1447 | ||
1448 | TVector3 gpos; | |
1449 | Float_t g[3]; | |
1450 | ||
1451 | // calculate new cluster position | |
596697e6 | 1452 | |
b43bdd14 | 1453 | recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr); |
6fdd30c4 | 1454 | recPoint->GetGlobalPosition(gpos); |
1455 | gpos.GetXYZ(g); | |
1456 | ||
1457 | // create a new cluster | |
596697e6 | 1458 | |
b43bdd14 | 1459 | (*fCaloClusterArr)[j] = new AliAODCaloCluster() ; |
1460 | AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ; | |
6a8be5c3 | 1461 | j++; |
6fdd30c4 | 1462 | clus->SetType(AliVCluster::kEMCALClusterv1); |
6544055e | 1463 | clus->SetE(clusterE); |
6fdd30c4 | 1464 | clus->SetPosition(g); |
f6e45fe0 | 1465 | clus->SetNCells(ncellsTrue); |
6fdd30c4 | 1466 | clus->SetCellsAbsId(absIds); |
1467 | clus->SetCellsAmplitudeFraction(ratios); | |
6fdd30c4 | 1468 | clus->SetChi2(-1); //not yet implemented |
bd9e8ebd | 1469 | clus->SetTOF(recPoint->GetTime()) ; //time-of-flight |
6fdd30c4 | 1470 | clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima |
5994e71f | 1471 | clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower()); |
b43bdd14 | 1472 | |
f46af216 | 1473 | if(ncells == ncellsTrue) |
1474 | { | |
f053064e | 1475 | Float_t elipAxis[2]; |
1476 | recPoint->GetElipsAxis(elipAxis); | |
1477 | clus->SetM02(elipAxis[0]*elipAxis[0]) ; | |
1478 | clus->SetM20(elipAxis[1]*elipAxis[1]) ; | |
1479 | clus->SetDispersion(recPoint->GetDispersion()); | |
1480 | } | |
f46af216 | 1481 | else if(fSelectCell) |
1482 | { | |
afaaef51 | 1483 | // In case some cells rejected, in unfolding case, recalculate |
1484 | // shower shape parameters and position | |
d20bc070 | 1485 | if(DebugLevel() > 1) |
f46af216 | 1486 | printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue); |
b43bdd14 | 1487 | |
f053064e | 1488 | AliVCaloCells* cells = 0x0; |
1489 | if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells(); | |
1490 | else cells = InputEvent()->GetEMCALCells(); | |
b43bdd14 | 1491 | |
f053064e | 1492 | fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus); |
1493 | fRecoUtils->RecalculateClusterPID(clus); | |
afaaef51 | 1494 | fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus); |
f053064e | 1495 | |
f053064e | 1496 | } |
6060ed91 | 1497 | |
6edb4448 | 1498 | // MC |
1499 | ||
1500 | if ( fSetCellMCLabelFromCluster == 1 ) SetClustersMCLabelFrom2SelectedLabels(recPoint,clus) ; | |
1501 | else if( fSetCellMCLabelFromCluster == 2 ) SetClustersMCLabelFromOriginalClusters(clus) ; | |
1502 | else | |
1503 | { | |
1504 | // Normal case, trust what the clusterizer has found | |
1505 | Int_t parentMult = 0; | |
1506 | Int_t *parentList = recPoint->GetParents(parentMult); | |
1507 | clus->SetLabel(parentList, parentMult); | |
1b5f58b6 | 1508 | // printf("Label list : "); |
1509 | // for(Int_t ilabel = 0; ilabel < parentMult; ilabel++ ) printf(" %d ",parentList[ilabel]); | |
1510 | // printf("\n"); | |
6edb4448 | 1511 | } |
1512 | ||
1513 | } // recPoints loop | |
1514 | ||
1515 | } | |
1516 | ||
1517 | //___________________________________________________________________________ | |
1518 | void AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs(Int_t & label) | |
1519 | { | |
1520 | // MC label for Cells not remapped after ESD filtering, do it here. | |
e59a6a2e | 1521 | |
6edb4448 | 1522 | if(label < 0) return ; |
1523 | ||
1524 | AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ; | |
1525 | if(!evt) return ; | |
1526 | ||
1527 | TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ; | |
1528 | if(!arr) return ; | |
6060ed91 | 1529 | |
6edb4448 | 1530 | if(label < arr->GetEntriesFast()) |
1531 | { | |
1532 | AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label)); | |
1533 | if(!particle) return ; | |
1534 | ||
1535 | if(label == particle->Label()) return ; // label already OK | |
1536 | //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label()); | |
1537 | } | |
1538 | //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast()); | |
1539 | ||
1540 | // loop on the particles list and check if there is one with the same label | |
1541 | for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ ) | |
1542 | { | |
1543 | AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind)); | |
48d9e9e7 | 1544 | if(!particle) continue ; |
1545 | ||
6edb4448 | 1546 | if(label == particle->Label()) |
f46af216 | 1547 | { |
6edb4448 | 1548 | label = ind; |
1549 | //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label); | |
1550 | return; | |
1551 | } | |
1552 | } | |
1553 | ||
1554 | label = -1; | |
1555 | ||
1556 | //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n"); | |
e59a6a2e | 1557 | |
6edb4448 | 1558 | } |
1559 | ||
1560 | //________________________________________________ | |
1561 | void AliAnalysisTaskEMCALClusterize::ResetArrays() | |
1562 | { | |
1563 | // Reset arrays containing information for all possible cells | |
1564 | for(Int_t j = 0; j < 12672; j++) | |
1565 | { | |
1566 | fOrgClusterCellId[j] = -1; | |
1567 | fCellLabels[j] = -1; | |
1568 | fCellSecondLabels[j] = -1; | |
1569 | fCellTime[j] = 0.; | |
1570 | fCellMatchdEta[j] = -999; | |
1571 | fCellMatchdPhi[j] = -999; | |
1572 | } | |
1573 | } | |
1574 | ||
1575 | //_____________________________________________________________________________________________________ | |
1576 | void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint * recPoint, | |
1577 | AliAODCaloCluster * clus) | |
1578 | { | |
1579 | // Set the cluster MC label, the digizer was filled with most likely MC label for all cells in original cluster | |
1580 | // Now check the second most likely MC label and add it to the new cluster | |
1581 | ||
1582 | Int_t parentMult = 0; | |
1583 | Int_t *parentList = recPoint->GetParents(parentMult); | |
1584 | clus->SetLabel(parentList, parentMult); | |
1585 | ||
1586 | //Write the second major contributor to each MC cluster. | |
1587 | Int_t iNewLabel ; | |
1588 | for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ) | |
1589 | { | |
1590 | ||
1591 | Int_t idCell = clus->GetCellAbsId(iLoopCell) ; | |
1592 | if(idCell>=0) | |
1593 | { | |
1594 | iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label. | |
1595 | for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ ) | |
b43bdd14 | 1596 | { |
6edb4448 | 1597 | if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label. |
6ffa213f | 1598 | if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0; |
6edb4448 | 1599 | } |
1600 | if (iNewLabel == 1) | |
1601 | { | |
1602 | Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ; | |
1603 | for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ) | |
6ffa213f | 1604 | { |
6edb4448 | 1605 | newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ; |
6ffa213f | 1606 | } |
6edb4448 | 1607 | |
1608 | newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ; | |
1609 | clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ; | |
1610 | delete [] newLabelArray; | |
1611 | } | |
1612 | }//positive cell number | |
1613 | } | |
1614 | } | |
1615 | ||
1616 | //___________________________________________________________________________________________________ | |
1617 | void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster * clus) | |
1618 | { | |
1619 | // Get the original clusters that contribute to the new cluster, assign the labels of such clusters | |
1620 | // to the new cluster. | |
1621 | // Only approximatedly valid when output are V1 clusters, handle with care | |
e59a6a2e | 1622 | |
11e7733f | 1623 | TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ... |
6edb4448 | 1624 | clArray.Reset(); |
1625 | Int_t nClu = 0; | |
1626 | Int_t nLabTotOrg = 0; | |
4d0b0c70 | 1627 | Float_t emax = -1; |
1628 | Int_t idMax = -1; | |
11e7733f | 1629 | |
1630 | AliVEvent * event = fEvent; | |
1631 | ||
1632 | //In case of embedding the MC signal is in the added event | |
1633 | AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); | |
1634 | if(aodIH && aodIH->GetEventToMerge()) //Embedding | |
1635 | event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal | |
1636 | ||
1637 | ||
6edb4448 | 1638 | //Find the clusters that originally had the cells |
1639 | for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ) | |
1640 | { | |
1641 | Int_t idCell = clus->GetCellAbsId(iLoopCell) ; | |
1642 | ||
1643 | if(idCell>=0) | |
1644 | { | |
1645 | Int_t idCluster = fOrgClusterCellId[idCell]; | |
1646 | ||
1647 | Bool_t set = kTRUE; | |
1648 | for(Int_t icl =0; icl < nClu; icl++) | |
1649 | { | |
11e7733f | 1650 | if(((Int_t)clArray.GetAt(icl))==-1 ) continue; |
6edb4448 | 1651 | if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE; |
1652 | // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n", | |
1653 | // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set); | |
1654 | } | |
1655 | if( set && idCluster >= 0) | |
1656 | { | |
1657 | clArray.SetAt(idCluster,nClu++); | |
1658 | //printf("******** idCluster %d \n",idCluster); | |
11e7733f | 1659 | nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels(); |
4d0b0c70 | 1660 | |
1661 | //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster); | |
1662 | ||
1663 | //Search highest E cluster | |
11e7733f | 1664 | AliVCluster * clOrg = event->GetCaloCluster(idCluster); |
4d0b0c70 | 1665 | //printf("\t E %f\n",clOrg->E()); |
1666 | if(emax < clOrg->E()) | |
1667 | { | |
1668 | emax = clOrg->E(); | |
1669 | idMax = idCluster; | |
1670 | } | |
6edb4448 | 1671 | } |
37d2296c | 1672 | } |
6edb4448 | 1673 | }// cell loop |
11e7733f | 1674 | |
6edb4448 | 1675 | if(nClu==0 || nLabTotOrg == 0) |
1676 | { | |
11e7733f | 1677 | //if(clus->E() > 0.25) printf("AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters() - Check: N org clusters %d, n tot labels %d, cluster E %f, n cells %d\n",nClu,nLabTotOrg,clus->E(), clus->GetNCells()); |
6edb4448 | 1678 | //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]); |
1679 | //printf("\n"); | |
1680 | } | |
1681 | ||
4d0b0c70 | 1682 | // Put the first in the list the cluster with highest energy |
1683 | if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position | |
1684 | { | |
1685 | Int_t maxIndex = -1; | |
1686 | Int_t firstCluster = ((Int_t)clArray.GetAt(0)); | |
1687 | for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ ) | |
1688 | { | |
1689 | if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster; | |
1690 | } | |
1691 | ||
11e7733f | 1692 | if(firstCluster >=0 && idMax >=0) |
1693 | { | |
1694 | clArray.SetAt(idMax,0); | |
1695 | clArray.SetAt(firstCluster,maxIndex); | |
1696 | } | |
4d0b0c70 | 1697 | } |
11e7733f | 1698 | |
6edb4448 | 1699 | // Get the labels list in the original clusters, assign all to the new cluster |
1700 | TArrayI clMCArray(nLabTotOrg) ; | |
1701 | clMCArray.Reset(); | |
1702 | ||
1703 | Int_t nLabTot = 0; | |
1704 | for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ ) | |
1705 | { | |
f88771a3 | 1706 | Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster); |
4d0b0c70 | 1707 | //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster); |
11e7733f | 1708 | AliVCluster * clOrg = event->GetCaloCluster(idCluster); |
6edb4448 | 1709 | Int_t nLab = clOrg->GetNLabels(); |
37d2296c | 1710 | |
6edb4448 | 1711 | for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ ) |
1712 | { | |
1713 | Int_t lab = clOrg->GetLabelAt(iLab) ; | |
1714 | if(lab>=0) | |
1715 | { | |
1716 | Bool_t set = kTRUE; | |
1717 | //printf("\t \t Set Label %d \n", lab); | |
1718 | for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++) | |
1719 | { | |
1720 | if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE; | |
1721 | //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n", | |
1722 | // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set); | |
1723 | } | |
1724 | if( set ) clMCArray.SetAt(lab,nLabTot++); | |
1725 | } | |
1726 | } | |
1727 | }// cluster loop | |
1728 | ||
1729 | // Set the final list of labels | |
1730 | ||
1731 | clus->SetLabel(clMCArray.GetArray(), nLabTot); | |
1732 | ||
e59a6a2e | 1733 | // printf("Final list of labels for new cluster : \n"); |
1734 | // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++) | |
1735 | // { | |
1736 | // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]); | |
1737 | // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]); | |
1738 | // printf(" org %d ",label); | |
1739 | // RemapMCLabelForAODs(label); | |
1740 | // printf(" new %d \n",label); | |
1741 | // } | |
1742 | // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc)); | |
6544055e | 1743 | } |
1744 | ||
3769e0cb | 1745 | |
6544055e | 1746 | //____________________________________________________________ |
1747 | void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects() | |
1748 | { | |
1749 | // Init geometry, create list of output clusters | |
1750 | ||
6a797202 | 1751 | |
1752 | fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0); | |
1753 | ||
1754 | if(fOutputAODBranchName.Length()==0) | |
3769e0cb | 1755 | { |
6a797202 | 1756 | fOutputAODBranchName = "newEMCALClustersArray"; |
1757 | printf("Cluster branch name not set, set it to newEMCALClustersArray \n"); | |
6544055e | 1758 | } |
6a797202 | 1759 | |
1760 | fOutputAODBranch->SetName(fOutputAODBranchName); | |
1761 | ||
1762 | if( fFillAODFile ) | |
3769e0cb | 1763 | { |
6a797202 | 1764 | //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0); |
1765 | ||
1766 | ||
1767 | //fOutputAODBranch->SetOwner(kFALSE); | |
1768 | ||
1769 | AddAODBranch("TClonesArray", &fOutputAODBranch); | |
6544055e | 1770 | } |
b43bdd14 | 1771 | |
6fdd30c4 | 1772 | } |
6544055e | 1773 | |
1774 | //_______________________________________________________ | |
1775 | void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *) | |
1776 | { | |
b43bdd14 | 1777 | // Do clusterization event by event, execute different steps |
596697e6 | 1778 | // 1) Do some checks on the kind of events (ESD, AOD) or if some filtering is needed, initializations |
1779 | // load things and clear arrays | |
b43bdd14 | 1780 | // 2) Clusterize a) just unfolding existing clusters (fJustUnfold) |
1781 | // b) recluster cells | |
1782 | // + convert cells into digits (calibrating them) | |
1783 | // + recluster digits into recPoints with the chosen clusterizer (V1, V1+Unfold,V2, NxN) | |
1784 | // with methods in AliEMCALClusterizer | |
1785 | // + transform recPoints into CaloClusters | |
1786 | // 3) Do some final calculations in the found clusters (track-matching) and put them in an AOD branch | |
1787 | ||
596697e6 | 1788 | //------- |
b43bdd14 | 1789 | // Step 1 |
1790 | ||
b0f874d4 | 1791 | |
b43bdd14 | 1792 | //Remove the contents of AOD branch output list set in the previous event |
6544055e | 1793 | fOutputAODBranch->Clear("C"); |
596697e6 | 1794 | |
cb68426a | 1795 | LoadBranches(); |
1796 | ||
1797 | //Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested | |
1798 | //If we need a centrality bin, we select only those events in the corresponding bin. | |
1799 | if( GetCentrality() && fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 ) | |
1800 | { | |
8fa7608a | 1801 | Float_t cen = GetEventCentrality(); |
cb68426a | 1802 | if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin. |
1803 | } | |
b0f874d4 | 1804 | |
596697e6 | 1805 | // intermediate array with new clusters : init the array only once or clear from previous event |
1806 | if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000); | |
1807 | else fCaloClusterArr->Delete();//Clear("C"); it leaks? | |
1808 | ||
5aa88029 | 1809 | InitGeometry(); // only once, must be done before OADB, geo OADB accessed here |
1810 | ||
6544055e | 1811 | //Get the event, do some checks and settings |
1812 | CheckAndGetEvent() ; | |
1813 | ||
a1eeb25d | 1814 | if (!fEvent) |
1815 | { | |
6544055e | 1816 | if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry()); |
1817 | return ; | |
1818 | } | |
6a797202 | 1819 | |
3769e0cb | 1820 | //Init pointers, geometry, clusterizer, ocdb, aodb |
6544055e | 1821 | |
6544055e | 1822 | |
3769e0cb | 1823 | if(fAccessOCDB) AccessOCDB(); |
1824 | if(fAccessOADB) AccessOADB(); // only once | |
b43bdd14 | 1825 | |
3769e0cb | 1826 | InitClusterization(); |
6544055e | 1827 | |
596697e6 | 1828 | //------- |
b43bdd14 | 1829 | // Step 2 |
1830 | ||
6544055e | 1831 | // Make clusters |
1832 | if (fJustUnfold) ClusterUnfolding(); | |
1833 | else ClusterizeCells() ; | |
6544055e | 1834 | |
596697e6 | 1835 | //------- |
b43bdd14 | 1836 | // Step 3 |
6544055e | 1837 | |
afc22931 | 1838 | FillCaloClusterInEvent(); |
b65291bf | 1839 | |
6544055e | 1840 | } |
1841 |