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