AliPhysics  32b88a8 (32b88a8)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEMCALClusterizeFast.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 // --- Root ---
17 #include <TClonesArray.h>
18 #include <TGeoManager.h>
19 #include <TObjArray.h>
20 #include <TString.h>
21 #include <TTree.h>
22 #include <TArrayI.h>
23 
24 // --- AliRoot ---
25 #include "AliAODCaloCluster.h"
26 #include "AliAODEvent.h"
27 #include "AliAnalysisManager.h"
28 #include "AliCDBEntry.h"
29 #include "AliCDBManager.h"
30 #include "AliCaloCalibPedestal.h"
31 #include "AliEMCALAfterBurnerUF.h"
32 #include "AliEMCALCalibData.h"
33 #include "AliEMCALClusterizerNxN.h"
34 #include "AliEMCALClusterizerv1.h"
35 #include "AliEMCALClusterizerv2.h"
36 #include "AliEMCALClusterizerFixedWindow.h"
37 #include "AliEMCALDigit.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliEMCALRecParam.h"
40 #include "AliEMCALRecPoint.h"
41 #include "AliEMCALRecoUtils.h"
42 #include "AliESDEvent.h"
43 #include "AliInputEventHandler.h"
44 #include "AliLog.h"
45 
47 
49 
50 //________________________________________________________________________
53  fRun(-1),
54  fDigitsArr(0),
55  fClusterArr(0),
56  fRecParam(new AliEMCALRecParam),
57  fClusterizer(0),
58  fUnfolder(0),
59  fJustUnfold(kFALSE),
60  fGeomName(),
61  fGeomMatrixSet(kFALSE),
62  fLoadGeomMatrices(kFALSE),
63  fOCDBpath(),
64  fCalibData(0),
65  fPedestalData(0),
66  fOutputAODBranch(0),
67  fOutputAODBrName(),
68  fRecoUtils(0),
69  fLoadCalib(kFALSE),
70  fLoadPed(kFALSE),
71  fAttachClusters(kTRUE),
72  fSubBackground(kFALSE),
73  fNPhi(4),
74  fNEta(4),
75  fShiftPhi(2),
76  fShiftEta(2),
77  fTRUShift(0),
78  fInputCellType(kFEEData),
79  fTrackName(),
80  fCaloCellsName(),
81  fCaloClustersName("newCaloClusters"),
82  fDoUpdateCells(kFALSE),
83  fDoClusterize(kTRUE),
84  fClusterBadChannelCheck(kFALSE),
85  fRejectExoticClusters(kFALSE),
86  fRejectExoticCells(kFALSE),
87  fFiducial(kFALSE),
88  fDoNonLinearity(kFALSE),
89  fRecalDistToBadChannels(kFALSE),
90  fSetCellMCLabelFromCluster(0),
91  fSetCellMCLabelFromEdepFrac(0),
92  fCaloCells(0),
93  fCaloClusters(0),
94  fEsd(0),
95  fAod(0),
96  fGeom(0)
97 {
98  // Constructor
99 
100  for(Int_t i = 0; i < AliEMCALGeoParams::fgkEMCALModules; i++) fGeomMatrix[i] = 0 ;
101  for(Int_t j = 0; j < fgkTotalCellNumber; j++)
102  { fOrgClusterCellId[j] =-1; fCellLabels[j] =-1 ; }
103 }
104 
105 //________________________________________________________________________
107  AliAnalysisTaskSE(name),
108  fRun(-1),
109  fDigitsArr(0),
110  fClusterArr(0),
111  fRecParam(new AliEMCALRecParam),
112  fClusterizer(0),
113  fUnfolder(0),
114  fJustUnfold(kFALSE),
115  fGeomName(),
116  fGeomMatrixSet(kFALSE),
117  fLoadGeomMatrices(kFALSE),
118  fOCDBpath(),
119  fCalibData(0),
120  fPedestalData(0),
121  fOutputAODBranch(0),
122  fOutputAODBrName(),
123  fRecoUtils(0),
124  fLoadCalib(kFALSE),
125  fLoadPed(kFALSE),
126  fAttachClusters(kTRUE),
127  fSubBackground(kFALSE),
128  fNPhi(4),
129  fNEta(4),
130  fShiftPhi(2),
131  fShiftEta(2),
132  fTRUShift(0),
133  fInputCellType(kFEEData),
134  fTrackName(),
135  fCaloCellsName(),
136  fCaloClustersName("newCaloClusters"),
137  fDoUpdateCells(kFALSE),
138  fDoClusterize(kTRUE),
139  fClusterBadChannelCheck(kFALSE),
140  fRejectExoticClusters(kFALSE),
141  fRejectExoticCells(kFALSE),
142  fFiducial(kFALSE),
143  fDoNonLinearity(kFALSE),
144  fRecalDistToBadChannels(kFALSE),
145  fSetCellMCLabelFromCluster(0),
146  fSetCellMCLabelFromEdepFrac(0),
147  fCaloCells(0),
148  fCaloClusters(0),
149  fEsd(0),
150  fAod(0),
151  fGeom(0)
152 {
153  // Constructor
154 
155  fBranchNames = "ESD:AliESDHeader.,AliESDRun.,EMCALCells.,EMCALTrigger. AOD:header,emcalCells";
156 
157  for(Int_t i = 0; i < AliEMCALGeoParams::fgkEMCALModules; i++) fGeomMatrix[i] = 0 ;
158  for(Int_t j = 0; j < fgkTotalCellNumber; j++)
159  { fOrgClusterCellId[j] =-1; fCellLabels[j] =-1 ; }
160 }
161 
162 //________________________________________________________________________
164 {
165  // Destructor.
166 
167  delete fClusterizer;
168  delete fUnfolder;
169  delete fRecoUtils;
170  delete fRecParam;
171 }
172 
173 //________________________________________________________________________
175 {
176  // Create output objects.
177 
178  if (!fOutputAODBrName.IsNull()) {
179  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
181  AddAODBranch("TClonesArray", &fOutputAODBranch);
182  AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
183  }
184 }
185 
186 //________________________________________________________________________
188 {
189  // Main loop, called for each event
190 
191  // if no operation is requested, return
192  if (!fDoClusterize && !fDoUpdateCells)
193  return;
194 
195  // remove the contents of output list set in the previous event
196  if (fOutputAODBranch)
197  fOutputAODBranch->Clear("C");
198 
199  fEsd = dynamic_cast<AliESDEvent*>(InputEvent());
200  fAod = dynamic_cast<AliAODEvent*>(InputEvent());
201 
202  if (!fEsd&&!fAod) {
203  Error("UserExec","Event not available");
204  return;
205  }
206 
207  LoadBranches();
208 
209  UInt_t offtrigger = 0;
210  if (fEsd) {
211  UInt_t mask1 = fEsd->GetESDRun()->GetDetectorsInDAQ();
212  UInt_t mask2 = fEsd->GetESDRun()->GetDetectorsInReco();
213  Bool_t desc1 = (mask1 >> 18) & 0x1;
214  Bool_t desc2 = (mask2 >> 18) & 0x1;
215  if (desc1==0 || desc2==0) { //AliDAQ::OfflineModuleName(180=="EMCAL"
216  AliError(Form("EMCAL not in DAQ/RECO: %u (%u)/%u (%u)",
217  mask1, fEsd->GetESDRun()->GetDetectorsInReco(),
218  mask2, fEsd->GetESDRun()->GetDetectorsInDAQ()));
219  return;
220  }
221  AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
222  offtrigger = ((AliInputEventHandler*)(am->GetInputEventHandler()))->IsEventSelected();
223  } else {
224  offtrigger = ((AliVAODHeader*)fAod->GetHeader())->GetOfflineTrigger();
225  }
226 
227  if (!MCEvent()) {
228  if (offtrigger & AliVEvent::kFastOnly) {
229  AliError(Form("EMCAL not in fast only partition"));
230  return;
231  }
232  }
233 
234  Init();
235 
236  if (fJustUnfold) {
237  AliWarning("Unfolding not implemented");
238  return;
239  }
240 
241  FillDigitsArray();
242 
243  if (fDoClusterize)
244  Clusterize();
245 
246  if (fDoUpdateCells)
247  UpdateCells();
248 
250  return;
251 
252  UpdateClusters();
254 
257 
258 }
259 
260 //________________________________________________________________________
261 void AliAnalysisTaskEMCALClusterizeFast::CopyClusters(TClonesArray *orig, TClonesArray *dest)
262 {
263  const Int_t Ncls = orig->GetEntries();
264 
265  for(Int_t i=0; i < Ncls; ++i) {
266  AliVCluster *oc = static_cast<AliVCluster*>(orig->At(i));
267 
268  if (!oc)
269  continue;
270 
271  if (!oc->IsEMCAL())
272  continue;
273 
274  AliVCluster *dc = static_cast<AliVCluster*>(dest->New(i));
275  dc->SetType(AliVCluster::kEMCALClusterv1);
276  dc->SetE(oc->E());
277  Float_t pos[3] = {0};
278  oc->GetPosition(pos);
279  dc->SetPosition(pos);
280  dc->SetNCells(oc->GetNCells());
281  dc->SetCellsAbsId(oc->GetCellsAbsId());
282  dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction());
283  dc->SetID(oc->GetID());
284  dc->SetDispersion(oc->GetDispersion());
285  dc->SetEmcCpvDistance(-1);
286  dc->SetChi2(-1);
287  dc->SetTOF(oc->GetTOF()); //time-of-flight
288  dc->SetNExMax(oc->GetNExMax()); //number of local maxima
289  dc->SetM02(oc->GetM02());
290  dc->SetM20(oc->GetM20());
291  dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel());
292  dc->SetMCEnergyFraction(oc->GetMCEnergyFraction());
293 
294  //MC
295  UInt_t nlabels = oc->GetNLabels();
296  Int_t *labels = oc->GetLabels();
297 
298  if (nlabels == 0 || !labels)
299  continue;
300 
301  AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
302  if (esdClus) {
303  TArrayI parents(nlabels, labels);
304  esdClus->AddLabels(parents);
305  }
306  else {
307  AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
308  if (aodClus)
309  aodClus->SetLabel(labels, nlabels);
310  }
311  }
312 }
313 
314 //________________________________________________________________________
316 {
317  // Clusterize
318 
319  if (fSubBackground) {
320  fClusterizer->SetInputCalibrated(kTRUE);
321  fClusterizer->SetCalibrationParameters(0);
322  }
323 
324  fClusterizer->Digits2Clusters("");
325 
326  if (fSubBackground) {
327  if (fCalibData) {
328  fClusterizer->SetInputCalibrated(kFALSE);
329  fClusterizer->SetCalibrationParameters(fCalibData);
330  }
331  }
332 }
333 
334 //________________________________________________________________________
336 {
337  // Fill digits array
338 
339  fDigitsArr->Clear("C");
340  switch (fInputCellType) {
341 
342  case kFEEData :
343  case kFEEDataMCOnly :
344  case kFEEDataExcludeMC :
345  {
346  // In case of MC productions done before aliroot tag v5-02-Rev09
347  // assing the cluster label to all the cells belonging to this cluster
348  // very rough
349  // Copied and simplified from AliEMCALTenderSupply
351  {
352  for (Int_t i = 0; i < fgkTotalCellNumber; i++)
353  {
354  fCellLabels [i] =-1 ;
355  fOrgClusterCellId[i] =-1 ;
356  }
357 
358  Int_t nClusters = InputEvent()->GetNumberOfCaloClusters();
359  for (Int_t i = 0; i < nClusters; i++)
360  {
361  AliVCluster *clus = InputEvent()->GetCaloCluster(i);
362 
363  if (!clus) continue;
364 
365  if (!clus->IsEMCAL()) continue ;
366 
367  Int_t label = clus->GetLabel();
368  UShort_t * index = clus->GetCellsAbsId() ;
369 
370  for(Int_t icell=0; icell < clus->GetNCells(); icell++)
371  {
373  fCellLabels[index[icell]] = label;
374 
375  fOrgClusterCellId[index[icell]] = i ; // index of the original cluster
376  } // cell in cluster loop
377  } // cluster loop
378  }
379 
380  Double_t avgE = 0; // for background subtraction
381  const Int_t ncells = fCaloCells->GetNumberOfCells();
382  for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell)
383  {
384  Double_t cellAmplitude=0, cellTime=0, cellEFrac = 0;
385  Short_t cellNumber=0;
386  Int_t cellMCLabel=-1;
387  if (fCaloCells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
388  break;
389 
390  if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
391 
392  if (cellMCLabel > 0 && cellEFrac < 1e-6)
393  cellEFrac = 1;
394 
395  if (cellAmplitude < 1e-6 || cellNumber < 0)
396  continue;
397 
399  if (cellMCLabel <= 0)
400  continue;
401  else {
402  cellAmplitude *= cellEFrac;
403  cellEFrac = 1;
404  }
405  }
406  else if (fInputCellType == kFEEDataExcludeMC) {
407  if (cellMCLabel > 0)
408  continue;
409  else {
410  cellAmplitude *= 1 - cellEFrac;
411  cellEFrac = 0;
412  }
413  }
414 
415  if(!AcceptCell(cellNumber)) continue;
416 
417  AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
418  (Float_t)cellAmplitude, (Float_t)cellTime,
419  AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
420 
422  {
423  Float_t energy = cellAmplitude;
424  Float_t time = cellTime;
425  fClusterizer->Calibrate(energy,time,cellNumber);
426  digit->SetAmplitude(energy);
427  avgE += energy;
428  }
429 
430  // New way to set the cell MC labels,
431  // valid only for MC productions with aliroot > v5-07-21
432  if(fSetCellMCLabelFromEdepFrac && fOrgClusterCellId[cellNumber] >= 0) // index can be negative if noisy cell that did not form cluster
433  {
434  fCellLabels[cellNumber] = idigit;
435 
436  AliVCluster *clus = 0;
437  Int_t iclus = fOrgClusterCellId[cellNumber];
438 
439  if(iclus < 0)
440  {
441  AliInfo("Negative original cluster index, skip \n");
442  continue;
443  }
444 
445  clus = InputEvent()->GetCaloCluster(iclus);
446 
447  for(Int_t icluscell=0; icluscell < clus->GetNCells(); icluscell++ )
448  {
449  if(cellNumber != clus->GetCellAbsId(icluscell)) continue ;
450 
451  // Get the energy deposition fraction.
452  Float_t eDepFrac[4];
453  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
454 
455  // Select the MC label contributing, only if enough energy deposition
456  TArrayI labeArr(0);
457  TArrayF eDepArr(0);
458  Int_t nLabels = 0;
459  for(Int_t imc = 0; imc < 4; imc++)
460  {
461  if(eDepFrac[imc] > 0 && clus->GetNLabels() > imc)
462  {
463  nLabels++;
464 
465  labeArr.Set(nLabels);
466  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
467 
468  eDepArr.Set(nLabels);
469  eDepArr.AddAt(eDepFrac[imc]*cellAmplitude, nLabels-1);
470  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
471  }
472  }
473 
474  if(nLabels > 0)
475  {
476  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
477  }
478  }
479  }
480 
481 
482  idigit++;
483  }
484 
485  if (fSubBackground) {
486  avgE /= fGeom->GetNumberOfSuperModules()*48*24;
487  Int_t ndigis = fDigitsArr->GetEntries();
488  for (Int_t i = 0; i < ndigis; ++i) {
489  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
490  Double_t energy = digit->GetAmplitude() - avgE;
491  if (energy<=0.001) {
492  digit->SetAmplitude(0);
493  } else {
494  digit->SetAmplitude(energy);
495  }
496  }
497  }
498  }
499  break;
500 
501  case kPattern :
502  {
503  // Fill digits from a pattern
504  Int_t maxd = fGeom->GetNCells() / 4;
505  for (Int_t idigit = 0; idigit < maxd; idigit++){
506  if (idigit % 24 == 12) idigit += 12;
507  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
508  digit->SetId(idigit * 4);
509  digit->SetTime(600);
510  digit->SetTimeR(600);
511  digit->SetIndexInList(idigit);
512  digit->SetType(AliEMCALDigit::kHG);
513  digit->SetAmplitude(0.1);
514  }
515  }
516  break;
517 
518  case kL0FastORs :
519  case kL0FastORsTC :
520  case kL1FastORs :
521  {
522  // Fill digits from FastORs
523 
524  AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
525 
526  if (!triggers || !(triggers->GetEntries() > 0))
527  return;
528 
529  Int_t idigit = 0;
530  triggers->Reset();
531 
532  while ((triggers->Next())) {
533  Float_t L0Amplitude = 0;
534  triggers->GetAmplitude(L0Amplitude);
535 
536  if (L0Amplitude <= 0 && fInputCellType != kL1FastORs)
537  continue;
538 
539  Int_t L1Amplitude = 0;
540  triggers->GetL1TimeSum(L1Amplitude);
541 
542  if (L1Amplitude <= 0 && fInputCellType == kL1FastORs)
543  continue;
544 
545  Int_t triggerTime = 0;
546  Int_t ntimes = 0;
547  triggers->GetNL0Times(ntimes);
548 
549  if (ntimes < 1 && fInputCellType == kL0FastORsTC)
550  continue;
551 
552  if (ntimes > 0) {
553  Int_t trgtimes[25];
554  triggers->GetL0Times(trgtimes);
555  triggerTime = trgtimes[0];
556  }
557 
558  Int_t triggerCol = 0, triggerRow = 0;
559  triggers->GetPosition(triggerCol, triggerRow);
560 
561  Int_t find = -1;
562  fGeom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
563 
564  if (find < 0)
565  continue;
566 
567  Int_t cidx[4] = {-1};
568  Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
569 
570  if (!ret)
571  continue;
572 
573  Float_t triggerAmplitude = 0;
574 
575  if (fInputCellType == kL1FastORs) {
576  triggerAmplitude = 0.25 * L1Amplitude; // it will add 4 cells for 1 amplitude
577  }
578  else {
579  triggerAmplitude = L0Amplitude; // 10 bit truncated, so it is already divided by 4
580  }
581 
582  for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
583  Int_t triggerNumber = cidx[idxpos];
584  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
585  digit->SetId(triggerNumber);
586  digit->SetTime(triggerTime);
587  digit->SetTimeR(triggerTime);
588  digit->SetIndexInList(idigit);
589  digit->SetType(AliEMCALDigit::kHG);
590  digit->SetAmplitude(triggerAmplitude);
591  idigit++;
592  }
593  }
594  }
595  break;
596  }
597 }
598 
599 //________________________________________________________________________________________
601 
602  Bool_t accept = kTRUE;
603  if(fRejectExoticCells) {
604  //Remove exotic cells before making digits
605  Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
606  fRecoUtils->SwitchOnRejectExoticCell();//switch on and off
607  Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
608  Bool_t isEx = fRecoUtils->IsExoticCell(cellNumber, fCaloCells, bunchCrossNo);
609  if(isEx) accept = kFALSE;
610  if(!exRemoval) fRecoUtils->SwitchOffRejectExoticCell();//switch on and off
611  }
612  return accept;
613 }
614 
615 //________________________________________________________________________________________
617 {
618  // Go through clusters one by one and process separate correction
619  // as those were defined or not
620 
621  Int_t nclusters = fCaloClusters->GetEntriesFast();
622  for (Int_t icluster=0; icluster < nclusters; ++icluster) {
623  AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
624  if (!clust)
625  continue;
626  if (!clust->IsEMCAL())
627  continue;
628 
629  // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
631  // careful, the the ClusterContainsBadChannel is dependent on
632  // SwitchOnBadChannelsRemoval, switching it ON automatically
633  // and returning to original value after processing
634  Bool_t badRemoval = fRecoUtils->IsBadChannelsRemovalSwitchedOn();
635  fRecoUtils->SwitchOnBadChannelsRemoval();
636 
637  Bool_t badResult = fRecoUtils->ClusterContainsBadChannel(fGeom, clust->GetCellsAbsId(), clust->GetNCells());
638 
639  // switch the bad channels removal back
640  if (!badRemoval)
641  fRecoUtils->SwitchOffBadChannelsRemoval();
642 
643  if (badResult) {
644  delete fCaloClusters->RemoveAt(icluster);
645  continue; //TODO is it really needed to remove it? Or should we flag it?
646  }
647  }
648 
649  // REMOVE EXOTIC CLUSTERS -------------------------------------
650  // does process local cell recalibration energy and time without replacing
651  // the global cell values, in case of no cell recalib done yet
652  if (fRejectExoticClusters) {
653  // careful, the IsExoticCluster is dependent on
654  // SwitchOnRejectExoticCell, switching it ON automatically
655  // and returning to original value after processing
656  Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
657  fRecoUtils->SwitchOnRejectExoticCell();
658 
659  // get bunch crossing
660  Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
661 
662  Bool_t exResult = fRecoUtils->IsExoticCluster(clust, fCaloCells, bunchCrossNo);
663 
664  // switch the exotic channels removal back
665  if (!exRemoval)
666  fRecoUtils->SwitchOffRejectExoticCell();
667 
668  if (exResult) {
669  delete fCaloClusters->RemoveAt(icluster);
670  continue; //TODO is it really needed to remove it? Or should we flag it?
671  }
672  }
673 
674  // FIDUCIAL CUT -----------------------------------------------
675  if (fFiducial) {
676  // depends on SetNumberOfCellsFromEMCALBorder
677  // SwitchOnNoFiducialBorderInEMCALEta0
678  if (!fRecoUtils->CheckCellFiducialRegion(fGeom, clust, fCaloCells)){
679  delete fCaloClusters->RemoveAt(icluster);
680  continue; //TODO it would be nice to store the distance
681  }
682  }
683 
684  // NONLINEARITY -----------------------------------------------
685  if (fDoNonLinearity) {
686  Float_t correctedEnergy = fRecoUtils->CorrectClusterEnergyLinearity(clust);
687  clust->SetE(correctedEnergy);
688  }
689 
690  // DISTANCE TO BAD CHANNELS -----------------------------------
692  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, clust);
693  }
694 
695  fCaloClusters->Compress();
696 }
697 
698 //________________________________________________________________________________________
699 void AliAnalysisTaskEMCALClusterizeFast::TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
700 {
701  Float_t g[3]={0};
702  c->GetPosition(g);
703  TVector3 gpos(g);
704 
705  Double_t dEtaMin = 1e9;
706  Double_t dPhiMin = 1e9;
707  Int_t imin = -1;
708  Double_t ceta = gpos.Eta();
709  Double_t cphi = gpos.Phi();
710  const Int_t ntrks = tarr->GetEntries();
711  for(Int_t t = 0; t<ntrks; ++t) {
712  AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
713  if (!track)
714  continue;
715  const AliExternalTrackParam *outp = track->GetOuterParam();
716  if (!outp)
717  continue;
718  Double_t trkPos[3] = {0.,0.,0.};
719  if (!outp->GetXYZ(trkPos))
720  continue;
721  TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
722  Double_t veta = vec.Eta();
723  Double_t vphi = vec.Phi();
724  if(vphi<0)
725  vphi += 2*TMath::Pi();
726  if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
727  continue;
728  Double_t dR = vec.DeltaR(gpos);
729  if(dR > 25)
730  continue;
731  Float_t tmpEta=0, tmpPhi=0;
732  if (0) {
733  AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
734  Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
735  if (!ret)
736  continue;
737  } else {
738  tmpEta = ceta - veta;
739  tmpPhi = cphi - vphi;
740  }
741  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
742  dEtaMin = tmpEta;
743  dPhiMin = tmpPhi;
744  imin = t;
745  }
746  }
747  c->SetEmcCpvDistance(imin);
748  c->SetTrackDistance(dPhiMin, dEtaMin);
749 }
750 
751 //________________________________________________________________________________________
753 {
754  // Cluster energy, global position, cells and their amplitude fractions are restored.
755 
756  // tracks array for track/cluster matching
757  TClonesArray *tarr = 0;
758  if (!fTrackName.IsNull()) {
759  tarr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackName));
760  if (!tarr) {
761  AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
762  }
763  }
764 
765  const Int_t Ncls = fClusterArr->GetEntries();
766  AliDebug(1, Form("total no of clusters %d", Ncls));
767 
768  for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i)
769  {
770  AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
771 
772  Int_t ncellsTrue = 0;
773  const Int_t ncells = recpoint->GetMultiplicity();
774  UShort_t absIds[ncells];
775  Double32_t ratios[ncells];
776  Int_t *dlist = recpoint->GetDigitsList();
777  Float_t *elist = recpoint->GetEnergiesList();
778  Double_t mcEnergy = 0;
779 
780  for (Int_t c = 0; c < ncells; ++c)
781  {
782  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
783  absIds[ncellsTrue] = digit->GetId();
784  ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
785 
786  if (digit->GetIparent(1) > 0)
787  mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
788 
789  ++ncellsTrue;
790  }
791 
792  if (ncellsTrue < 1)
793  {
794  AliWarning("Skipping cluster with no cells");
795  continue;
796  }
797 
798  // calculate new cluster position
799  TVector3 gpos;
800  recpoint->GetGlobalPosition(gpos);
801  Float_t g[3];
802  gpos.GetXYZ(g);
803 
804  AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
805 
806  AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
807  c->SetType(AliVCluster::kEMCALClusterv1);
808  c->SetE(recpoint->GetEnergy());
809  c->SetPosition(g);
810  c->SetNCells(ncellsTrue);
811  c->SetCellsAbsId(absIds);
812  c->SetCellsAmplitudeFraction(ratios);
813  c->SetID(recpoint->GetUniqueID());
814  c->SetDispersion(recpoint->GetDispersion());
815  c->SetEmcCpvDistance(-1);
816  c->SetChi2(-1);
817  c->SetTOF(recpoint->GetTime()) ; //time-of-flight
818  c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
819  Float_t elipAxis[2];
820  recpoint->GetElipsAxis(elipAxis);
821  c->SetM02(elipAxis[0]*elipAxis[0]);
822  c->SetM20(elipAxis[1]*elipAxis[1]);
823  c->SetMCEnergyFraction(mcEnergy);
824 
825  //
826  // MC labels
827  //
828  Int_t parentMult = 0;
829  Int_t *parentList = recpoint->GetParents(parentMult);
830  Float_t *parentListDE = recpoint->GetParentsDE(); // deposited energy
831 
832  c->SetLabel(parentList, parentMult);
833  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
834 
835  //
836  // Set the cell energy deposition fraction map:
837  //
838  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
839  {
840  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
841 
842  // Get the digit that originated this cell cluster
843  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
844 
845  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
846  {
847  Int_t idigit = fCellLabels[absIds[icell]];
848 
849  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
850 
851  // Find the 4 MC labels that contributed to the cluster and their
852  // deposited energy in the current digit
853 
854  mcEdepFracPerCell[icell] = 0; // init
855 
856  Int_t nparents = dig->GetNiparent();
857  if ( nparents > 0 )
858  {
859  Int_t digLabel =-1 ;
860  Float_t edep = 0 ;
861  Float_t edepTot = 0 ;
862  Float_t mcEDepFrac[4] = {0,0,0,0};
863 
864  // all parents in digit
865  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
866  {
867  digLabel = dig->GetIparent (jndex+1);
868  edep = dig->GetDEParent(jndex+1);
869  edepTot += edep;
870 
871  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
872  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
873  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
874  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
875  } // all prarents in digit
876 
877  // Divide energy deposit by total deposited energy
878  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
879  if(edepTot > 0.01)
880  {
881  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
882  }
883  } // at least one parent label in digit
884  } // cell in cluster loop
885 
886  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
887 
888  delete [] mcEdepFracPerCell;
889 
890  } // at least one parent in cluster, do the cell primary packing
891 
892  //
893  // Track matching
894  //
895  if (tarr)
896  TrackClusterMatching(c, tarr);
897  }
898 }
899 
900 //________________________________________________________________________
902 {
903  // Update cells in case re-calibration was done.
904  if (!fCalibData&&!fSubBackground)
905  return;
906 
907  const Int_t ncells = fCaloCells->GetNumberOfCells();
908  const Int_t ndigis = fDigitsArr->GetEntries();
909  if (ncells!=ndigis) {
910  fCaloCells->DeleteContainer();
911  fCaloCells->CreateContainer(ndigis);
912  }
913  for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
914  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
915  Double_t cellAmplitude = digit->GetCalibAmp();
916  Short_t cellNumber = digit->GetId();
917  Double_t cellTime = digit->GetTime();
918  fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
919  }
920 }
921 
922 //________________________________________________________________________
924 {
925  // Update cells in case re-calibration was done.
926 
927  const Int_t nents = fCaloClusters->GetEntries();
928  for (Int_t i=0;i<nents;++i) {
929  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
930  if (!c)
931  continue;
932  if (c->IsEMCAL())
933  delete fCaloClusters->RemoveAt(i);
934  }
935 
936  fCaloClusters->Compress();
937 
939 }
940 
941 //________________________________________________________________________________________
943 {
944  // Select clusterization/unfolding algorithm and set all the needed parameters.
945 
946  if (InputEvent()->GetRunNumber()==fRun)
947  return;
948  fRun = InputEvent()->GetRunNumber();
949 
950  if (fJustUnfold){
951  // init the unfolding afterburner
952  delete fUnfolder;
953  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
954  return;
955  }
956 
957  if (fGeomName.Length()>0)
958  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
959  else
960  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
961  if (!fGeom) {
962  AliFatal("Geometry not available!!!");
963  return;
964  }
965 
966  if (!fGeomMatrixSet) {
967  if (fLoadGeomMatrices) {
968  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
969  if (fGeomMatrix[mod]){
970  if (DebugLevel() > 2)
971  fGeomMatrix[mod]->Print();
972  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
973  }
974  }
975  } else { // get matrix from file (work around bug in aliroot)
976  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
977  const TGeoHMatrix *gm = 0;
978  if (fEsd) {
979  gm = fEsd->GetEMCALMatrix(mod);
980  } else {
981  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
982  if(!aodheader) AliFatal("Not a standard AOD");
983  if (aodheader) {
984  gm = aodheader->GetEMCALMatrix(mod);
985  }
986  }
987  if (gm) {
988  if (DebugLevel() > 2)
989  gm->Print();
990  fGeom->SetMisalMatrix(gm,mod);
991  }
992  }
993  }
994  fGeomMatrixSet=kTRUE;
995  }
996 
997  // setup digit array if needed
998  if (!fDigitsArr) {
999  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
1000  fDigitsArr->SetOwner(1);
1001  }
1002 
1003  // then setup clusterizer
1004  if (fClusterizer) {
1005  // avoid to delete digits array
1006  fClusterizer->SetDigitsArr(0);
1007  delete fClusterizer;
1008  }
1009  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1010  fClusterizer = new AliEMCALClusterizerv1(fGeom);
1011  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
1012  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
1013  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1014  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1015  fClusterizer = clusterizer;
1016  }
1017  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1018  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1019  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
1020  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
1021  clusterizer->SetNphi(fNPhi);
1022  clusterizer->SetNeta(fNEta);
1023  clusterizer->SetShiftPhi(fShiftPhi);
1024  clusterizer->SetShiftEta(fShiftEta);
1025  clusterizer->SetTRUshift(fTRUShift);
1026  fClusterizer = clusterizer;
1027  }
1028  else {
1029  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1030  }
1031  fClusterizer->InitParameters(fRecParam);
1032 
1033  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
1034  AliCDBManager *cdb = AliCDBManager::Instance();
1035  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
1036  cdb->SetDefaultStorage(fOCDBpath);
1037  if (fRun!=cdb->GetRun())
1038  cdb->SetRun(fRun);
1039  }
1040  if (!fCalibData&&fLoadCalib&&fRun>0) {
1041  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
1042  if (entry)
1043  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
1044  if (!fCalibData)
1045  AliFatal("Calibration parameters not found in CDB!");
1046  }
1047  if (!fPedestalData&&fLoadPed&&fRun>0) {
1048  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
1049  if (entry)
1050  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
1051  }
1052  if (fCalibData) {
1053  fClusterizer->SetInputCalibrated(kFALSE);
1054  fClusterizer->SetCalibrationParameters(fCalibData);
1055  } else {
1056  fClusterizer->SetInputCalibrated(kTRUE);
1057  }
1058  fClusterizer->SetCaloCalibPedestal(fPedestalData);
1059  fClusterizer->SetJustClusters(kTRUE);
1060  fClusterizer->SetDigitsArr(fDigitsArr);
1061  fClusterizer->SetOutput(0);
1062  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1063 
1064  // Get the emcal cells
1066  if (fCaloCellsName.IsNull()) {
1067  fCaloCells = InputEvent()->GetEMCALCells();
1068  }
1069  else {
1070  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
1071  if (!fCaloCells)
1072  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
1073  }
1074  if (!fCaloCells)
1075  AliFatal("Could not get EMCal cells!");
1076  }
1077 
1078  // Set output clusters collection
1079  if (!fAttachClusters) {
1081  return;
1082  }
1083 
1084  if (!fCaloClusters) {
1085  if (fCaloClustersName.IsNull()) { //overwrite mode
1086  if (fEsd)
1087  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
1088  else if (fAod)
1089  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
1090  }
1091  else {
1092  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
1093 
1094  if (!fCaloClusters) {
1095  if (fEsd)
1096  fCaloClusters = new TClonesArray("AliESDCaloCluster");
1097  else if (fAod)
1098  fCaloClusters = new TClonesArray("AliAODCaloCluster");
1099 
1100  fCaloClusters->SetName(fCaloClustersName);
1101  InputEvent()->AddObject(fCaloClusters);
1102  }
1103  }
1104 
1105  if (!fCaloClusters)
1106  AliFatal("Could not get cluster collection!");
1107 
1108  TClass *cl = fCaloClusters->GetClass();
1109  if (!cl->GetBaseClass("AliVCluster")) {
1110  AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
1111  fCaloClusters = 0;
1112  return;
1113  }
1114  }
1115 }
void Print(std::ostream &o, const char *name, Double_t dT, Double_t dVM, Double_t alldT, Double_t alldVM)
Definition: PlotSysInfo.C:121
double Double_t
Definition: External.C:58
TString fOutputAODBrName
AOD Branch with output clusters.
TCanvas * c
Definition: TestFitELoss.C:172
AliEMCALAfterBurnerUF * fUnfolder
clusterizer
virtual void TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
virtual void CopyClusters(TClonesArray *orig, TClonesArray *dest)
short Short_t
Definition: External.C:23
energy
TGeoHMatrix * fGeomMatrix[AliEMCALGeoParams::fgkEMCALModules]
ClassImp(AliAnalysisTaskEMCALClusterizeFast) AliAnalysisTaskEMCALClusterizeFast
unsigned short UShort_t
Definition: External.C:28
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53