AliPhysics  d497547 (d497547)
 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 
827  c->SetClusterMCEdepFractionFromEdepArray(parentListDE);
828 
829  //
830  // Set the cell energy deposition fraction map:
831  //
832  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
833  {
834  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
835 
836  // Get the digit that originated this cell cluster
837  //AliVCaloCells* cells = InputEvent()->GetEMCALCells();
838 
839  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
840  {
841  Int_t idigit = fCellLabels[absIds[icell]];
842 
843  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
844 
845  // Find the 4 MC labels that contributed to the cluster and their
846  // deposited energy in the current digit
847 
848  mcEdepFracPerCell[icell] = 0; // init
849 
850  Int_t nparents = dig->GetNiparent();
851  if ( nparents > 0 )
852  {
853  Int_t digLabel =-1 ;
854  Float_t edep = 0 ;
855  Float_t edepTot = 0 ;
856  Float_t mcEDepFrac[4] = {0,0,0,0};
857 
858  // all parents in digit
859  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
860  {
861  digLabel = dig->GetIparent (jndex+1);
862  edep = dig->GetDEParent(jndex+1);
863  edepTot += edep;
864 
865  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
866  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
867  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
868  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
869  } // all prarents in digit
870 
871  // Divide energy deposit by total deposited energy
872  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
873  if(edepTot > 0.01)
874  {
875  mcEdepFracPerCell[icell] = c->PackMCEdepFraction(mcEDepFrac);
876  }
877  } // at least one parent label in digit
878  } // cell in cluster loop
879 
880  c->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
881 
882  delete [] mcEdepFracPerCell;
883 
884  } // at least one parent in cluster, do the cell primary packing
885 
886  //
887  // Track matching
888  //
889  if (tarr)
890  TrackClusterMatching(c, tarr);
891  }
892 }
893 
894 //________________________________________________________________________
896 {
897  // Update cells in case re-calibration was done.
898  if (!fCalibData&&!fSubBackground)
899  return;
900 
901  const Int_t ncells = fCaloCells->GetNumberOfCells();
902  const Int_t ndigis = fDigitsArr->GetEntries();
903  if (ncells!=ndigis) {
904  fCaloCells->DeleteContainer();
905  fCaloCells->CreateContainer(ndigis);
906  }
907  for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
908  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
909  Double_t cellAmplitude = digit->GetCalibAmp();
910  Short_t cellNumber = digit->GetId();
911  Double_t cellTime = digit->GetTime();
912  fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
913  }
914 }
915 
916 //________________________________________________________________________
918 {
919  // Update cells in case re-calibration was done.
920 
921  const Int_t nents = fCaloClusters->GetEntries();
922  for (Int_t i=0;i<nents;++i) {
923  AliVCluster *c = static_cast<AliVCluster*>(fCaloClusters->At(i));
924  if (!c)
925  continue;
926  if (c->IsEMCAL())
927  delete fCaloClusters->RemoveAt(i);
928  }
929 
930  fCaloClusters->Compress();
931 
933 }
934 
935 //________________________________________________________________________________________
937 {
938  // Select clusterization/unfolding algorithm and set all the needed parameters.
939 
940  if (InputEvent()->GetRunNumber()==fRun)
941  return;
942  fRun = InputEvent()->GetRunNumber();
943 
944  if (fJustUnfold){
945  // init the unfolding afterburner
946  delete fUnfolder;
947  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
948  return;
949  }
950 
951  if (fGeomName.Length()>0)
952  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
953  else
954  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
955  if (!fGeom) {
956  AliFatal("Geometry not available!!!");
957  return;
958  }
959 
960  if (!fGeomMatrixSet) {
961  if (fLoadGeomMatrices) {
962  for(Int_t mod=0; mod < fGeom->GetNumberOfSuperModules(); ++mod) {
963  if (fGeomMatrix[mod]){
964  if (DebugLevel() > 2)
965  fGeomMatrix[mod]->Print();
966  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod);
967  }
968  }
969  } else { // get matrix from file (work around bug in aliroot)
970  for(Int_t mod=0; mod < fGeom->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
971  const TGeoHMatrix *gm = 0;
972  if (fEsd) {
973  gm = fEsd->GetEMCALMatrix(mod);
974  } else {
975  AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(fAod->GetHeader());
976  if(!aodheader) AliFatal("Not a standard AOD");
977  if (aodheader) {
978  gm = aodheader->GetEMCALMatrix(mod);
979  }
980  }
981  if (gm) {
982  if (DebugLevel() > 2)
983  gm->Print();
984  fGeom->SetMisalMatrix(gm,mod);
985  }
986  }
987  }
988  fGeomMatrixSet=kTRUE;
989  }
990 
991  // setup digit array if needed
992  if (!fDigitsArr) {
993  fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
994  fDigitsArr->SetOwner(1);
995  }
996 
997  // then setup clusterizer
998  if (fClusterizer) {
999  // avoid to delete digits array
1000  fClusterizer->SetDigitsArr(0);
1001  delete fClusterizer;
1002  }
1003  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1004  fClusterizer = new AliEMCALClusterizerv1(fGeom);
1005  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN) {
1006  AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom);
1007  clusterizer->SetNRowDiff(fRecParam->GetNRowDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1008  clusterizer->SetNColDiff(fRecParam->GetNColDiff()); //MV: already done in AliEMCALClusterizer::InitParameters
1009  fClusterizer = clusterizer;
1010  }
1011  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1012  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1013  else if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerFW) {
1014  AliEMCALClusterizerFixedWindow *clusterizer = new AliEMCALClusterizerFixedWindow(fGeom);
1015  clusterizer->SetNphi(fNPhi);
1016  clusterizer->SetNeta(fNEta);
1017  clusterizer->SetShiftPhi(fShiftPhi);
1018  clusterizer->SetShiftEta(fShiftEta);
1019  clusterizer->SetTRUshift(fTRUShift);
1020  fClusterizer = clusterizer;
1021  }
1022  else {
1023  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1024  }
1025  fClusterizer->InitParameters(fRecParam);
1026 
1027  if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
1028  AliCDBManager *cdb = AliCDBManager::Instance();
1029  if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
1030  cdb->SetDefaultStorage(fOCDBpath);
1031  if (fRun!=cdb->GetRun())
1032  cdb->SetRun(fRun);
1033  }
1034  if (!fCalibData&&fLoadCalib&&fRun>0) {
1035  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
1036  if (entry)
1037  fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
1038  if (!fCalibData)
1039  AliFatal("Calibration parameters not found in CDB!");
1040  }
1041  if (!fPedestalData&&fLoadPed&&fRun>0) {
1042  AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
1043  if (entry)
1044  fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
1045  }
1046  if (fCalibData) {
1047  fClusterizer->SetInputCalibrated(kFALSE);
1048  fClusterizer->SetCalibrationParameters(fCalibData);
1049  } else {
1050  fClusterizer->SetInputCalibrated(kTRUE);
1051  }
1052  fClusterizer->SetCaloCalibPedestal(fPedestalData);
1053  fClusterizer->SetJustClusters(kTRUE);
1054  fClusterizer->SetDigitsArr(fDigitsArr);
1055  fClusterizer->SetOutput(0);
1056  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1057 
1058  // Get the emcal cells
1060  if (fCaloCellsName.IsNull()) {
1061  fCaloCells = InputEvent()->GetEMCALCells();
1062  }
1063  else {
1064  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
1065  if (!fCaloCells)
1066  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
1067  }
1068  if (!fCaloCells)
1069  AliFatal("Could not get EMCal cells!");
1070  }
1071 
1072  // Set output clusters collection
1073  if (!fAttachClusters) {
1075  return;
1076  }
1077 
1078  if (!fCaloClusters) {
1079  if (fCaloClustersName.IsNull()) { //overwrite mode
1080  if (fEsd)
1081  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
1082  else if (fAod)
1083  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
1084  }
1085  else {
1086  fCaloClusters = static_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloClustersName));
1087 
1088  if (!fCaloClusters) {
1089  if (fEsd)
1090  fCaloClusters = new TClonesArray("AliESDCaloCluster");
1091  else if (fAod)
1092  fCaloClusters = new TClonesArray("AliAODCaloCluster");
1093 
1094  fCaloClusters->SetName(fCaloClustersName);
1095  InputEvent()->AddObject(fCaloClusters);
1096  }
1097  }
1098 
1099  if (!fCaloClusters)
1100  AliFatal("Could not get cluster collection!");
1101 
1102  TClass *cl = fCaloClusters->GetClass();
1103  if (!cl->GetBaseClass("AliVCluster")) {
1104  AliFatal(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloClusters->GetName()));
1105  fCaloClusters = 0;
1106  return;
1107  }
1108  }
1109 }
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.
energy
Definition: HFPtSpectrum.C:44
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
TGeoHMatrix * fGeomMatrix[AliEMCALGeoParams::fgkEMCALModules]
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