AliPhysics  1adf5bd (1adf5bd)
 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  Float_t cellAmplitude=0;
385  Double_t cellTime=0, amp = 0, cellEFrac = 0;
386  Short_t cellNumber=0;
387  Int_t cellMCLabel=-1;
388  if (fCaloCells->GetCell(icell, cellNumber, amp, cellTime, cellMCLabel, cellEFrac) != kTRUE)
389  break;
390 
391  cellAmplitude = amp; // compilation problem
392 
393  if (fSetCellMCLabelFromCluster) cellMCLabel = fCellLabels[cellNumber];
394 
395  if (cellMCLabel > 0 && cellEFrac < 1e-6)
396  cellEFrac = 1;
397 
398  if (cellAmplitude < 1e-6 || cellNumber < 0)
399  continue;
400 
402  if (cellMCLabel <= 0)
403  continue;
404  else {
405  cellAmplitude *= cellEFrac;
406  cellEFrac = 1;
407  }
408  }
409  else if (fInputCellType == kFEEDataExcludeMC) {
410  if (cellMCLabel > 0)
411  continue;
412  else {
413  cellAmplitude *= 1 - cellEFrac;
414  cellEFrac = 0;
415  }
416  }
417 
418  if(!AcceptCell(cellNumber)) continue;
419 
420  //
421  // New way to set the cell MC labels,
422  // valid only for MC productions with aliroot > v5-07-21
423  //
424  TArrayI labeArr(0);
425  TArrayF eDepArr(0);
426  Int_t nLabels = 0;
427 
428  if(fSetCellMCLabelFromEdepFrac && fOrgClusterCellId[cellNumber] >= 0) // index can be negative if noisy cell that did not form cluster
429  {
430  cellMCLabel = -1;
431 
432  fCellLabels[cellNumber] = idigit;
433 
434  Int_t iclus = fOrgClusterCellId[cellNumber];
435 
436  if(iclus < 0)
437  {
438  AliInfo("Negative original cluster index, skip \n");
439  continue;
440  }
441 
442  AliVCluster *clus = InputEvent()->GetCaloCluster(iclus);
443 
444  for(Int_t icluscell=0; icluscell < clus->GetNCells(); icluscell++ )
445  {
446  if(cellNumber != clus->GetCellAbsId(icluscell)) continue ;
447 
448  // Select the MC label contributing, only if enough energy deposition
449 
450  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(cellNumber, clus, MCEvent(),cellAmplitude, labeArr, eDepArr);
451  nLabels = labeArr.GetSize();
452  }
453  }
454 
455  AliEMCALDigit *digit = new((*fDigitsArr)[idigit]) AliEMCALDigit(cellMCLabel, cellMCLabel, cellNumber,
456  cellAmplitude, (Float_t)cellTime,
457  AliEMCALDigit::kHG,idigit, 0, 0, cellEFrac*cellAmplitude);
458 
459  if(nLabels > 0)
460  {
461  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
462  }
463 
464 
466  {
467  Float_t energy = cellAmplitude;
468  Float_t time = cellTime;
469  fClusterizer->Calibrate(energy,time,cellNumber);
470  digit->SetAmplitude(energy);
471  avgE += energy;
472  }
473 
474  idigit++;
475  }
476 
477  if (fSubBackground) {
478  avgE /= fGeom->GetNumberOfSuperModules()*48*24;
479  Int_t ndigis = fDigitsArr->GetEntries();
480  for (Int_t i = 0; i < ndigis; ++i) {
481  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
482  Double_t energy = digit->GetAmplitude() - avgE;
483  if (energy<=0.001) {
484  digit->SetAmplitude(0);
485  } else {
486  digit->SetAmplitude(energy);
487  }
488  }
489  }
490  }
491  break;
492 
493  case kPattern :
494  {
495  // Fill digits from a pattern
496  Int_t maxd = fGeom->GetNCells() / 4;
497  for (Int_t idigit = 0; idigit < maxd; idigit++){
498  if (idigit % 24 == 12) idigit += 12;
499  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
500  digit->SetId(idigit * 4);
501  digit->SetTime(600);
502  digit->SetTimeR(600);
503  digit->SetIndexInList(idigit);
504  digit->SetType(AliEMCALDigit::kHG);
505  digit->SetAmplitude(0.1);
506  }
507  }
508  break;
509 
510  case kL0FastORs :
511  case kL0FastORsTC :
512  case kL1FastORs :
513  {
514  // Fill digits from FastORs
515 
516  AliVCaloTrigger *triggers = InputEvent()->GetCaloTrigger("EMCAL");
517 
518  if (!triggers || !(triggers->GetEntries() > 0))
519  return;
520 
521  Int_t idigit = 0;
522  triggers->Reset();
523 
524  while ((triggers->Next())) {
525  Float_t L0Amplitude = 0;
526  triggers->GetAmplitude(L0Amplitude);
527 
528  if (L0Amplitude <= 0 && fInputCellType != kL1FastORs)
529  continue;
530 
531  Int_t L1Amplitude = 0;
532  triggers->GetL1TimeSum(L1Amplitude);
533 
534  if (L1Amplitude <= 0 && fInputCellType == kL1FastORs)
535  continue;
536 
537  Int_t triggerTime = 0;
538  Int_t ntimes = 0;
539  triggers->GetNL0Times(ntimes);
540 
541  if (ntimes < 1 && fInputCellType == kL0FastORsTC)
542  continue;
543 
544  if (ntimes > 0) {
545  Int_t trgtimes[25];
546  triggers->GetL0Times(trgtimes);
547  triggerTime = trgtimes[0];
548  }
549 
550  Int_t triggerCol = 0, triggerRow = 0;
551  triggers->GetPosition(triggerCol, triggerRow);
552 
553  Int_t find = -1;
554  fGeom->GetAbsFastORIndexFromPositionInEMCAL(triggerCol, triggerRow, find);
555 
556  if (find < 0)
557  continue;
558 
559  Int_t cidx[4] = {-1};
560  Bool_t ret = fGeom->GetCellIndexFromFastORIndex(find, cidx);
561 
562  if (!ret)
563  continue;
564 
565  Float_t triggerAmplitude = 0;
566 
567  if (fInputCellType == kL1FastORs) {
568  triggerAmplitude = 0.25 * L1Amplitude; // it will add 4 cells for 1 amplitude
569  }
570  else {
571  triggerAmplitude = L0Amplitude; // 10 bit truncated, so it is already divided by 4
572  }
573 
574  for (Int_t idxpos = 0; idxpos < 4; idxpos++) {
575  Int_t triggerNumber = cidx[idxpos];
576  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
577  digit->SetId(triggerNumber);
578  digit->SetTime(triggerTime);
579  digit->SetTimeR(triggerTime);
580  digit->SetIndexInList(idigit);
581  digit->SetType(AliEMCALDigit::kHG);
582  digit->SetAmplitude(triggerAmplitude);
583  idigit++;
584  }
585  }
586  }
587  break;
588  }
589 }
590 
591 //________________________________________________________________________________________
593 
594  Bool_t accept = kTRUE;
595  if(fRejectExoticCells) {
596  //Remove exotic cells before making digits
597  Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
598  fRecoUtils->SwitchOnRejectExoticCell();//switch on and off
599  Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
600  Bool_t isEx = fRecoUtils->IsExoticCell(cellNumber, fCaloCells, bunchCrossNo);
601  if(isEx) accept = kFALSE;
602  if(!exRemoval) fRecoUtils->SwitchOffRejectExoticCell();//switch on and off
603  }
604  return accept;
605 }
606 
607 //________________________________________________________________________________________
609 {
610  // Go through clusters one by one and process separate correction
611  // as those were defined or not
612 
613  Int_t nclusters = fCaloClusters->GetEntriesFast();
614  for (Int_t icluster=0; icluster < nclusters; ++icluster) {
615  AliVCluster *clust = static_cast<AliVCluster*>(fCaloClusters->At(icluster));
616  if (!clust)
617  continue;
618  if (!clust->IsEMCAL())
619  continue;
620 
621  // REMOVE CLUSTERS WITH BAD CELLS -----------------------------
623  // careful, the the ClusterContainsBadChannel is dependent on
624  // SwitchOnBadChannelsRemoval, switching it ON automatically
625  // and returning to original value after processing
626  Bool_t badRemoval = fRecoUtils->IsBadChannelsRemovalSwitchedOn();
627  fRecoUtils->SwitchOnBadChannelsRemoval();
628 
629  Bool_t badResult = fRecoUtils->ClusterContainsBadChannel(fGeom, clust->GetCellsAbsId(), clust->GetNCells());
630 
631  // switch the bad channels removal back
632  if (!badRemoval)
633  fRecoUtils->SwitchOffBadChannelsRemoval();
634 
635  if (badResult) {
636  delete fCaloClusters->RemoveAt(icluster);
637  continue; //TODO is it really needed to remove it? Or should we flag it?
638  }
639  }
640 
641  // REMOVE EXOTIC CLUSTERS -------------------------------------
642  // does process local cell recalibration energy and time without replacing
643  // the global cell values, in case of no cell recalib done yet
644  if (fRejectExoticClusters) {
645  // careful, the IsExoticCluster is dependent on
646  // SwitchOnRejectExoticCell, switching it ON automatically
647  // and returning to original value after processing
648  Bool_t exRemoval = fRecoUtils->IsRejectExoticCell();
649  fRecoUtils->SwitchOnRejectExoticCell();
650 
651  // get bunch crossing
652  Int_t bunchCrossNo = InputEvent()->GetBunchCrossNumber();
653 
654  Bool_t exResult = fRecoUtils->IsExoticCluster(clust, fCaloCells, bunchCrossNo);
655 
656  // switch the exotic channels removal back
657  if (!exRemoval)
658  fRecoUtils->SwitchOffRejectExoticCell();
659 
660  if (exResult) {
661  delete fCaloClusters->RemoveAt(icluster);
662  continue; //TODO is it really needed to remove it? Or should we flag it?
663  }
664  }
665 
666  // FIDUCIAL CUT -----------------------------------------------
667  if (fFiducial) {
668  // depends on SetNumberOfCellsFromEMCALBorder
669  // SwitchOnNoFiducialBorderInEMCALEta0
670  if (!fRecoUtils->CheckCellFiducialRegion(fGeom, clust, fCaloCells)){
671  delete fCaloClusters->RemoveAt(icluster);
672  continue; //TODO it would be nice to store the distance
673  }
674  }
675 
676  // NONLINEARITY -----------------------------------------------
677  if (fDoNonLinearity) {
678  Float_t correctedEnergy = fRecoUtils->CorrectClusterEnergyLinearity(clust);
679  clust->SetE(correctedEnergy);
680  }
681 
682  // DISTANCE TO BAD CHANNELS -----------------------------------
684  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, clust);
685  }
686 
687  fCaloClusters->Compress();
688 }
689 
690 //________________________________________________________________________________________
691 void AliAnalysisTaskEMCALClusterizeFast::TrackClusterMatching(AliVCluster *c, TClonesArray *tarr)
692 {
693  Float_t g[3]={0};
694  c->GetPosition(g);
695  TVector3 gpos(g);
696 
697  Double_t dEtaMin = 1e9;
698  Double_t dPhiMin = 1e9;
699  Int_t imin = -1;
700  Double_t ceta = gpos.Eta();
701  Double_t cphi = gpos.Phi();
702  const Int_t ntrks = tarr->GetEntries();
703  for(Int_t t = 0; t<ntrks; ++t) {
704  AliVTrack *track = static_cast<AliVTrack*>(tarr->At(t));
705  if (!track)
706  continue;
707  const AliExternalTrackParam *outp = track->GetOuterParam();
708  if (!outp)
709  continue;
710  Double_t trkPos[3] = {0.,0.,0.};
711  if (!outp->GetXYZ(trkPos))
712  continue;
713  TVector3 vec(trkPos[0],trkPos[1],trkPos[2]);
714  Double_t veta = vec.Eta();
715  Double_t vphi = vec.Phi();
716  if(vphi<0)
717  vphi += 2*TMath::Pi();
718  if (TMath::Abs(veta)>0.75 || (vphi<70*TMath::DegToRad()) || (vphi>190*TMath::DegToRad()))
719  continue;
720  Double_t dR = vec.DeltaR(gpos);
721  if(dR > 25)
722  continue;
723  Float_t tmpEta=0, tmpPhi=0;
724  if (0) {
725  AliExternalTrackParam trkParTemp(*outp); // retrieve the starting point every time before the extrapolation
726  Bool_t ret = fRecoUtils->ExtrapolateTrackToCluster(&trkParTemp, c, fRecoUtils->GetMass(), fRecoUtils->GetStep(), tmpEta, tmpPhi);
727  if (!ret)
728  continue;
729  } else {
730  tmpEta = ceta - veta;
731  tmpPhi = cphi - vphi;
732  }
733  if (TMath::Abs(tmpEta)<TMath::Abs(dEtaMin) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMin)) {
734  dEtaMin = tmpEta;
735  dPhiMin = tmpPhi;
736  imin = t;
737  }
738  }
739  c->SetEmcCpvDistance(imin);
740  c->SetTrackDistance(dPhiMin, dEtaMin);
741 }
742 
743 //________________________________________________________________________________________
745 {
746  // Cluster energy, global position, cells and their amplitude fractions are restored.
747 
748  // tracks array for track/cluster matching
749  TClonesArray *tarr = 0;
750  if (!fTrackName.IsNull()) {
751  tarr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackName));
752  if (!tarr) {
753  AliError(Form("Cannot get tracks named %s", fTrackName.Data()));
754  }
755  }
756 
757  const Int_t Ncls = fClusterArr->GetEntries();
758  AliDebug(1, Form("total no of clusters %d", Ncls));
759 
760  for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i)
761  {
762  AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
763 
764  Int_t ncellsTrue = 0;
765  const Int_t ncells = recpoint->GetMultiplicity();
766  UShort_t absIds[ncells];
767  Double32_t ratios[ncells];
768  Int_t *dlist = recpoint->GetDigitsList();
769  Float_t *elist = recpoint->GetEnergiesList();
770  Double_t mcEnergy = 0;
771 
772  for (Int_t c = 0; c < ncells; ++c)
773  {
774  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
775  absIds[ncellsTrue] = digit->GetId();
776  ratios[ncellsTrue] = elist[c]/digit->GetAmplitude();
777 
778  if (digit->GetIparent(1) > 0)
779  mcEnergy += digit->GetDEParent(1)/recpoint->GetEnergy();
780 
781  ++ncellsTrue;
782  }
783 
784  if (ncellsTrue < 1)
785  {
786  AliWarning("Skipping cluster with no cells");
787  continue;
788  }
789 
790  // calculate new cluster position
791  TVector3 gpos;
792  recpoint->GetGlobalPosition(gpos);
793  Float_t g[3];
794  gpos.GetXYZ(g);
795 
796  AliDebug(1, Form("energy %f", recpoint->GetEnergy()));
797 
798  AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
799  c->SetType(AliVCluster::kEMCALClusterv1);
800  c->SetE(recpoint->GetEnergy());
801  c->SetPosition(g);
802  c->SetNCells(ncellsTrue);
803  c->SetCellsAbsId(absIds);
804  c->SetCellsAmplitudeFraction(ratios);
805  c->SetID(recpoint->GetUniqueID());
806  c->SetDispersion(recpoint->GetDispersion());
807  c->SetEmcCpvDistance(-1);
808  c->SetChi2(-1);
809  c->SetTOF(recpoint->GetTime()) ; //time-of-flight
810  c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
811  Float_t elipAxis[2];
812  recpoint->GetElipsAxis(elipAxis);
813  c->SetM02(elipAxis[0]*elipAxis[0]);
814  c->SetM20(elipAxis[1]*elipAxis[1]);
815  c->SetMCEnergyFraction(mcEnergy);
816 
817  //
818  // MC labels
819  //
820  Int_t parentMult = 0;
821  Int_t *parentList = recpoint->GetParents(parentMult);
822  Float_t *parentListDE = recpoint->GetParentsDE(); // deposited energy
823 
824  c->SetLabel(parentList, parentMult);
825  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
826 
827  //
828  // Set the cell energy deposition fraction map:
829  //
830  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
831  {
832  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
833 
834  // Get the digit that originated this cell cluster
835  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
836 
837  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
838  {
839  Int_t idigit = fCellLabels[absIds[icell]];
840 
841  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
842 
843  // Find the 4 MC labels that contributed to the cluster and their
844  // deposited energy in the current digit
845 
846  mcEdepFracPerCell[icell] = 0; // init
847 
848  Int_t nparents = dig->GetNiparent();
849  if ( nparents > 0 )
850  {
851  Int_t digLabel =-1 ;
852  Float_t edep = 0 ;
853  Float_t edepTot = 0 ;
854  Float_t mcEDepFrac[4] = {0,0,0,0};
855 
856  // all parents in digit
857  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
858  {
859  digLabel = dig->GetIparent (jndex+1);
860  edep = dig->GetDEParent(jndex+1);
861  edepTot += edep;
862 
863  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
864  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
865  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
866  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
867  } // all prarents in digit
868 
869  // Divide energy deposit by total deposited energy
870  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
871  if(edepTot > 0.01)
872  {
873  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
874  }
875  } // at least one parent label in digit
876  } // cell in cluster loop
877 
878  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
879 
880  delete [] mcEdepFracPerCell;
881 
882  } // at least one parent in cluster, do the cell primary packing
883 
884  //
885  // Track matching
886  //
887  if (tarr)
888  TrackClusterMatching(c, tarr);
889  }
890 }
891 
892 //________________________________________________________________________
894 {
895  // Update cells in case re-calibration was done.
896  if (!fCalibData&&!fSubBackground)
897  return;
898 
899  const Int_t ncells = fCaloCells->GetNumberOfCells();
900  const Int_t ndigis = fDigitsArr->GetEntries();
901  if (ncells!=ndigis) {
902  fCaloCells->DeleteContainer();
903  fCaloCells->CreateContainer(ndigis);
904  }
905  for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
906  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
907  Double_t cellAmplitude = digit->GetCalibAmp();
908  Short_t cellNumber = digit->GetId();
909  Double_t cellTime = digit->GetTime();
910  fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
911  }
912 }
913 
914 //________________________________________________________________________
916 {
917  // Update cells in case re-calibration was done.
918 
919  const Int_t nents = fCaloClusters->GetEntries();
920  for (Int_t i=0;i<nents;++i) {
921  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
922  if (!c)
923  continue;
924  if (c->IsEMCAL())
925  delete fCaloClusters->RemoveAt(i);
926  }
927 
928  fCaloClusters->Compress();
929 
931 }
932 
933 //________________________________________________________________________________________
935 {
936  // Select clusterization/unfolding algorithm and set all the needed parameters.
937 
938  if (InputEvent()->GetRunNumber()==fRun)
939  return;
940  fRun = InputEvent()->GetRunNumber();
941 
942  if (fJustUnfold){
943  // init the unfolding afterburner
944  delete fUnfolder;
945  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
946  return;
947  }
948 
949  if (fGeomName.Length()>0)
950  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
951  else
952  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
953  if (!fGeom) {
954  AliFatal("Geometry not available!!!");
955  return;
956  }
957 
958  if (!fGeomMatrixSet) {
959  if (fLoadGeomMatrices) {
960  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
961  if (fGeomMatrix[mod]){
962  if (DebugLevel() > 2)
963  fGeomMatrix[mod]->Print();
964  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
965  }
966  }
967  } else { // get matrix from file (work around bug in aliroot)
968  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
969  const TGeoHMatrix *gm = 0;
970  if (fEsd) {
971  gm = fEsd->GetEMCALMatrix(mod);
972  } else {
973  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
974  if(!aodheader) AliFatal("Not a standard AOD");
975  if (aodheader) {
976  gm = aodheader->GetEMCALMatrix(mod);
977  }
978  }
979  if (gm) {
980  if (DebugLevel() > 2)
981  gm->Print();
982  fGeom->SetMisalMatrix(gm,mod);
983  }
984  }
985  }
986  fGeomMatrixSet=kTRUE;
987  }
988 
989  // setup digit array if needed
990  if (!fDigitsArr) {
991  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
992  fDigitsArr->SetOwner(1);
993  }
994 
995  // then setup clusterizer
996  if (fClusterizer) {
997  // avoid to delete digits array
998  fClusterizer->SetDigitsArr(0);
999  delete fClusterizer;
1000  }
1001  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1002  fClusterizer = new AliEMCALClusterizerv1(fGeom);
1003  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
1004  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
1005  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1006  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1007  fClusterizer = clusterizer;
1008  }
1009  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1010  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1011  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
1012  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
1013  clusterizer->SetNphi(fNPhi);
1014  clusterizer->SetNeta(fNEta);
1015  clusterizer->SetShiftPhi(fShiftPhi);
1016  clusterizer->SetShiftEta(fShiftEta);
1017  clusterizer->SetTRUshift(fTRUShift);
1018  fClusterizer = clusterizer;
1019  }
1020  else {
1021  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1022  }
1023  fClusterizer->InitParameters(fRecParam);
1024 
1025  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
1026  AliCDBManager *cdb = AliCDBManager::Instance();
1027  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
1028  cdb->SetDefaultStorage(fOCDBpath);
1029  if (fRun!=cdb->GetRun())
1030  cdb->SetRun(fRun);
1031  }
1032  if (!fCalibData&&fLoadCalib&&fRun>0) {
1033  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
1034  if (entry)
1035  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
1036  if (!fCalibData)
1037  AliFatal("Calibration parameters not found in CDB!");
1038  }
1039  if (!fPedestalData&&fLoadPed&&fRun>0) {
1040  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
1041  if (entry)
1042  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
1043  }
1044  if (fCalibData) {
1045  fClusterizer->SetInputCalibrated(kFALSE);
1046  fClusterizer->SetCalibrationParameters(fCalibData);
1047  } else {
1048  fClusterizer->SetInputCalibrated(kTRUE);
1049  }
1050  fClusterizer->SetCaloCalibPedestal(fPedestalData);
1051  fClusterizer->SetJustClusters(kTRUE);
1052  fClusterizer->SetDigitsArr(fDigitsArr);
1053  fClusterizer->SetOutput(0);
1054  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1055 
1056  // Get the emcal cells
1058  if (fCaloCellsName.IsNull()) {
1059  fCaloCells = InputEvent()->GetEMCALCells();
1060  }
1061  else {
1062  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
1063  if (!fCaloCells)
1064  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
1065  }
1066  if (!fCaloCells)
1067  AliFatal("Could not get EMCal cells!");
1068  }
1069 
1070  // Set output clusters collection
1071  if (!fAttachClusters) {
1073  return;
1074  }
1075 
1076  if (!fCaloClusters) {
1077  if (fCaloClustersName.IsNull()) { //overwrite mode
1078  if (fEsd)
1079  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
1080  else if (fAod)
1081  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
1082  }
1083  else {
1084  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
1085 
1086  if (!fCaloClusters) {
1087  if (fEsd)
1088  fCaloClusters = new TClonesArray("AliESDCaloCluster");
1089  else if (fAod)
1090  fCaloClusters = new TClonesArray("AliAODCaloCluster");
1091 
1092  fCaloClusters->SetName(fCaloClustersName);
1093  InputEvent()->AddObject(fCaloClusters);
1094  }
1095  }
1096 
1097  if (!fCaloClusters)
1098  AliFatal("Could not get cluster collection!");
1099 
1100  TClass *cl = fCaloClusters->GetClass();
1101  if (!cl->GetBaseClass("AliVCluster")) {
1102  AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
1103  fCaloClusters = 0;
1104  return;
1105  }
1106  }
1107 }
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