AliPhysics  8dc8609 (8dc8609)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEMCALClusterize.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 <TString.h>
18 #include <TRefArray.h>
19 #include <TClonesArray.h>
20 #include <TTree.h>
21 #include <TGeoManager.h>
22 #include <TROOT.h>
23 #include <TInterpreter.h>
24 #include <TFile.h>
25 
26 // --- AliRoot Analysis Steering
27 #include "AliAnalysisTask.h"
28 #include "AliAnalysisManager.h"
29 #include "AliESDEvent.h"
30 #include "AliGeomManager.h"
31 #include "AliVCaloCells.h"
32 #include "AliAODCaloCluster.h"
33 #include "AliCDBManager.h"
34 #include "AliCDBStorage.h"
35 #include "AliCDBEntry.h"
36 #include "AliLog.h"
37 #include "AliVEventHandler.h"
38 #include "AliAODInputHandler.h"
39 #include "AliOADBContainer.h"
40 #include "AliAODMCParticle.h"
41 #include "AliCentrality.h"
42 #include "AliMultSelection.h"
43 
44 // --- EMCAL
45 #include "AliEMCALAfterBurnerUF.h"
46 #include "AliEMCALGeometry.h"
47 #include "AliEMCALClusterizerNxN.h"
48 #include "AliEMCALClusterizerv1.h"
49 #include "AliEMCALClusterizerv2.h"
50 #include "AliEMCALRecPoint.h"
51 #include "AliEMCALDigit.h"
52 
53 #include "AliAnalysisTaskEMCALClusterize.h"
54 
56 ClassImp(AliAnalysisTaskEMCALClusterize) ;
58 
59 //______________________________________________________________________________
60 // Constructor.
61 //______________________________________________________________________________
62 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
63 : AliAnalysisTaskSE(name)
64 , fEvent(0)
65 , fGeom(0), fGeomName("")
66 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
67 , fOCDBpath(""), fAccessOCDB(kFALSE)
68 , fDigitsArr(0), fClusterArr(0)
69 , fCaloClusterArr(0), fCaloCells(0)
70 , fRecParam(0), fClusterizer(0)
71 , fUnfolder(0), fJustUnfold(kFALSE)
72 , fOutputAODBranch(0), fOutputAODBranchName("")
73 , fOutputAODCells (0), fOutputAODCellsName (""), fInputCaloCellsName ("")
74 , fOutputAODBranchSet(0)
75 , fFillAODFile(kFALSE), fFillAODHeader(0)
76 , fFillAODCaloCells(0), fRun(-1)
77 , fRecoUtils(0), fConfigName("")
78 , fOrgClusterCellId()
79 , fCellLabels(), fCellSecondLabels(), fCellTime()
80 , fCellMatchdEta(), fCellMatchdPhi()
81 , fRecalibrateWithClusterTime(0)
82 , fMaxEvent(0), fDoTrackMatching(kFALSE), fUpdateCell(0)
83 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
84 , fRejectBelowThreshold(kFALSE)
85 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
86 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
87 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
88 , fConstantTimeShift(0)
89 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
90 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
91 , fSetCellMCLabelFromCluster(0)
92 , fSetCellMCLabelFromEdepFrac(0)
93 , fRemapMCLabelForAODs(0)
94 , fInputFromFilter(0)
95 , fTCardCorrEmulation(0), fTCardCorrClusEnerConserv(0)
96 , fRandom(0), fRandomizeTCard(1)
97 , fTCardCorrMinAmp(0.01), fTCardCorrMaxInduced(100)
98 , fPrintOnce(0)
99 
100 {
101  for(Int_t i = 0; i < 22; i++)
102  {
103  fGeomMatrix[i] = 0;
104  fTCardCorrInduceEnerProb [i] = 0;
105  fTCardCorrInduceEnerFracMax[i] = 100;
106  fTCardCorrInduceEnerFracMin[i] =-100;
107 
108  for(Int_t j = 0; j < 4 ; j++)
109  {
110  fTCardCorrInduceEner [j][i] = 0 ;
111  fTCardCorrInduceEnerFrac [j][i] = 0 ;
112  fTCardCorrInduceEnerFracP1 [j][i] = 0 ;
113  fTCardCorrInduceEnerFracWidth[j][i] = 0 ;
114  }
115  }
116 
117  ResetArrays();
118 
119  fCentralityBin[0] = fCentralityBin[1]=-1;
120 }
121 
122 //______________________________________________________________
124 //______________________________________________________________
125 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
126 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
127 , fEvent(0)
128 , fGeom(0), fGeomName("")
129 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
130 , fOCDBpath(""), fAccessOCDB(kFALSE)
131 , fDigitsArr(0), fClusterArr(0)
132 , fCaloClusterArr(0), fCaloCells(0)
133 , fRecParam(0), fClusterizer(0)
134 , fUnfolder(0), fJustUnfold(kFALSE)
135 , fOutputAODBranch(0), fOutputAODBranchName("")
136 , fOutputAODCells (0), fOutputAODCellsName (""), fInputCaloCellsName ("")
137 , fOutputAODBranchSet(0)
138 , fFillAODFile(kFALSE), fFillAODHeader(0)
139 , fFillAODCaloCells(0), fRun(-1)
140 , fRecoUtils(0), fConfigName("")
141 , fOrgClusterCellId()
142 , fCellLabels(), fCellSecondLabels(), fCellTime()
143 , fCellMatchdEta(), fCellMatchdPhi()
144 , fRecalibrateWithClusterTime(0)
145 , fMaxEvent(0), fDoTrackMatching(kFALSE), fUpdateCell(0)
146 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
147 , fRejectBelowThreshold(kFALSE)
148 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
149 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
150 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
151 , fConstantTimeShift(0)
152 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
153 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
154 , fSetCellMCLabelFromCluster(0)
155 , fSetCellMCLabelFromEdepFrac(0)
156 , fRemapMCLabelForAODs(0)
157 , fInputFromFilter(0)
158 , fTCardCorrEmulation(0), fTCardCorrClusEnerConserv(0)
159 , fRandom(0), fRandomizeTCard(1)
160 , fTCardCorrMinAmp(0.01), fTCardCorrMaxInduced(100)
161 , fPrintOnce(0)
162 {
163  for(Int_t i = 0; i < 22; i++)
164  {
165  fGeomMatrix[i] = 0;
166  fTCardCorrInduceEnerProb [i] = 0;
167  fTCardCorrInduceEnerFracMax[i] = 100;
168  fTCardCorrInduceEnerFracMin[i] =-100;
169 
170  for(Int_t j = 0; j < 4 ; j++)
171  {
172  fTCardCorrInduceEner [j][i] = 0 ;
173  fTCardCorrInduceEnerFrac [j][i] = 0 ;
174  fTCardCorrInduceEnerFracP1 [j][i] = 0 ;
175  fTCardCorrInduceEnerFracWidth[j][i] = 0 ;
176  }
177  }
178 
179  ResetArrays();
180 
181  fCentralityBin[0] = fCentralityBin[1]=-1;
182 }
183 
184 //_______________________________________________________________
186 //_______________________________________________________________
187 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
188 {
189  if (fDigitsArr)
190  {
191  fDigitsArr->Clear("C");
192  delete fDigitsArr;
193  }
194 
195  if (fClusterArr)
196  {
197  fClusterArr->Delete();
198  delete fClusterArr;
199  }
200 
201  if (fCaloClusterArr)
202  {
203  fCaloClusterArr->Delete();
204  delete fCaloClusterArr;
205  }
206 
207  if(fClusterizer) delete fClusterizer;
208  if(fUnfolder) delete fUnfolder;
209  if(fRecoUtils) delete fRecoUtils;
210 }
211 
212 //_______________________________________________________
220 //_______________________________________________________________________________
221 Bool_t AliAnalysisTaskEMCALClusterize::AcceptCell(Int_t absID)
222 {
223  if ( absID < 0 || absID >= 24*48*fGeom->GetNumberOfSuperModules() )
224  return kFALSE;
225 
226  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
227  if (!fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
228  return kFALSE;
229 
230  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
231 
232  // Do not include bad channels found in analysis,
233  if ( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
234  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi) )
235  return kFALSE;
236 
237  return kTRUE;
238 }
239 
240 //_______________________________________________________
244 //_______________________________________________________
245 Bool_t AliAnalysisTaskEMCALClusterize::AcceptEventEMCAL()
246 {
247  if(!fSelectEMCALEvent) return kTRUE; // accept
248 
249  if(fEMCALEnergyCut <= 0) return kTRUE; // accept
250 
251  Int_t nCluster = fEvent -> GetNumberOfCaloClusters();
252  Int_t bc = fEvent -> GetBunchCrossNumber();
253 
254  for(Int_t icalo = 0; icalo < nCluster; icalo++)
255  {
256  AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
257 
258  if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
259  fRecoUtils->IsGoodCluster(clus,fGeom,fCaloCells,bc))
260  {
261 
262  AliDebug(1, Form("Accept : E %2.2f > %2.2f, nCells %d > %d",
263  clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut));
264 
265  return kTRUE;
266  }
267 
268  }// loop
269 
270  AliDebug(1,"Reject");
271 
272  return kFALSE;
273 }
274 
275 //_______________________________________________
278 //_______________________________________________
279 void AliAnalysisTaskEMCALClusterize::AccessOADB()
280 {
281  // Set it only once, unless run changed
282  if ( fOADBSet ) return ;
283 
284  TString pass = GetPass();
285 
286  AliInfo(Form("Get AODB parameters from EMCAL in %s for run %d, and <%s>",fOADBFilePath.Data(),fRun,pass.Data()));
287 
288  Int_t nSM = fGeom->GetNumberOfSuperModules();
289 
290  // Bad map
291  if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
292  {
293  AliOADBContainer *contBC=new AliOADBContainer("");
294  contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
295 
296  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRun);
297 
298  if(arrayBC)
299  {
300  fRecoUtils->SwitchOnDistToBadChannelRecalculation();
301  AliInfo("Remove EMCAL bad cells");
302 
303  for (Int_t i=0; i < nSM; ++i)
304  {
305  TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
306 
307  if (hbm)
308  delete hbm;
309 
310  hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
311 
312  if (!hbm)
313  {
314  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
315  continue;
316  }
317 
318  hbm->SetDirectory(0);
319  fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
320 
321  } // loop
322  } else AliInfo("Do NOT remove EMCAL bad channels"); // run array
323 
324  delete contBC;
325  } // Remove bad
326 
327  // Energy Recalibration
328  if(fRecoUtils->IsRecalibrationOn())
329  {
330  AliOADBContainer *contRF=new AliOADBContainer("");
331 
332  contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
333 
334  TObjArray *recal=(TObjArray*)contRF->GetObject(fRun);
335 
336  if(recal)
337  {
338  TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
339 
340  if(recalpass)
341  {
342  TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
343 
344  if(recalib)
345  {
346  AliInfo("Recalibrate EMCAL");
347  for (Int_t i=0; i<nSM; ++i)
348  {
349  TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
350 
351  if (h)
352  delete h;
353 
354  h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
355 
356  if (!h)
357  {
358  AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
359  continue;
360  }
361 
362  h->SetDirectory(0);
363 
364  fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
365  } // SM loop
366  } else AliInfo("Do NOT recalibrate EMCAL, no params object array"); // array ok
367  } else AliInfo("Do NOT recalibrate EMCAL, no params for pass"); // array pass ok
368  } else AliInfo("Do NOT recalibrate EMCAL, no params for run"); // run number array ok
369 
370  delete contRF;
371  } // Recalibration on
372 
373  // Energy Recalibration, apply on top of previous calibration factors
374  if(fRecoUtils->IsRunDepRecalibrationOn())
375  {
376  AliOADBContainer *contRFTD=new AliOADBContainer("");
377 
378  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
379 
380  TH1S *htd=(TH1S*)contRFTD->GetObject(fRun);
381 
382  //If it did not exist for this run, get closes one
383  if (!htd)
384  {
385  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",fRun));
386  // let's get the closest fRun instead then..
387  Int_t lower = 0;
388  Int_t ic = 0;
389  Int_t maxEntry = contRFTD->GetNumberOfEntries();
390 
391  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < fRun) ) {
392  lower = ic;
393  ic++;
394  }
395 
396  Int_t closest = lower;
397  if ( (ic<maxEntry) &&
398  (contRFTD->LowerLimit(ic)-fRun) < (fRun - contRFTD->UpperLimit(lower)) ) {
399  closest = ic;
400  }
401 
402  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
403  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
404  }
405 
406  if(htd)
407  {
408  AliInfo("Recalibrate (Temperature) EMCAL");
409 
410  for (Int_t ism=0; ism<nSM; ++ism)
411  {
412  for (Int_t icol=0; icol<48; ++icol)
413  {
414  for (Int_t irow=0; irow<24; ++irow)
415  {
416  Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
417 
418  Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
419  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
420  //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
421  // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
422  fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
423  } // columns
424  } // rows
425  } // SM loop
426  } else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
427 
428  delete contRFTD;
429  } // Run by Run T calibration
430 
431  // Time Recalibration
432  if(fRecoUtils->IsTimeRecalibrationOn())
433  {
434  AliOADBContainer *contTRF=new AliOADBContainer("");
435 
436  contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
437 
438  TObjArray *trecal=(TObjArray*)contTRF->GetObject(fRun);
439 
440  if(trecal)
441  {
442  // pass number should be pass1 except on Run1 and special cases
443  TString passM = pass;
444  if ( pass=="spc_calo" ) passM = "pass3";
445  if ( fRun > 209121 ) passM = "pass1"; // run2 periods
446  if ( pass == "muon_calo_pass1" && fRun > 209121 && fRun < 244284 )
447  passM = "pass0";//period LHC15a-m
448 
449  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
450 
451  if(trecalpass)
452  {
453  AliInfo("Time Recalibrate EMCAL");
454  for (Int_t ibc = 0; ibc < 4; ++ibc)
455  {
456  TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
457 
458  if (h)
459  delete h;
460 
461  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
462 
463  if (!h)
464  {
465  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
466  continue;
467  }
468 
469  h->SetDirectory(0);
470 
471  fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
472  } // bunch crossing loop
473  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
474  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
475 
476  delete contTRF;
477  } // Time recalibration on
478 
479  // L1 Phase Time Recalibration
480  if( fRecoUtils->IsL1PhaseInTimeRecalibrationOn() )
481  {
482  // init default maps first
483  if (!fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationArray())
484  fRecoUtils->InitEMCALL1PhaseInTimeRecalibration() ;
485 
486  AliOADBContainer *contBC = new AliOADBContainer("");
487 
488  TFile *timeFile=new TFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"read");
489  if (!timeFile || timeFile->IsZombie())
490  {
491  AliFatal(Form("EMCALTimeL1PhaseCalib.root was not found in the path provided: %s",fOADBFilePath.Data()));
492  return ;
493  }
494 
495  if (timeFile) delete timeFile;
496 
497  contBC->InitFromFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"AliEMCALTimeL1PhaseCalib");
498 
499  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(fRun);
500  if (!arrayBC)
501  {
502  AliError(Form("No external L1 phase in time calibration set for run number: %d", fRun));
503  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
504  }
505  else
506  {
507  // Only 1 L1 phase correction possible, except special cases
508  TString pass2 = "pass1";
509 
510  if ( pass=="muon_calo_pass1" && fRun > 209121 && fRun < 244284 )
511  pass2 = "pass0"; // period LHC15a-m
512 
513  TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(pass2);
514  if (!arrayBCpass)
515  {
516  AliError(Form("No external L1 phase in time calibration set for: %d -%s", fRun,pass2.Data()));
517  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
518  }
519  else AliInfo("Recalibrate L1 Phase time");
520 
521  if(arrayBCpass)
522  {
523  if ( DebugLevel()>0 ) arrayBCpass->Print();
524 
525  TH1C *h = fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationForAllSM();
526  if (h) delete h;
527 
528  h = (TH1C*)arrayBCpass->FindObject(Form("h%d",fRun));
529 
530  if (!h)
531  {
532  AliFatal(Form("There is no calibration histogram h%d for this run",fRun));
533  return;
534  }
535 
536  h->SetDirectory(0);
537  fRecoUtils->SetEMCALL1PhaseInTimeRecalibrationForAllSM(h);
538  }
539  }
540 
541  delete contBC;
542  } // L1 Phase Time Recalibration
543 
544  // Parameters already set once, so do not it again, unless run changes
545  fOADBSet = kTRUE;
546 }
547 
548 //_________________________________________________
551 //_________________________________________________
552 void AliAnalysisTaskEMCALClusterize::AccessOCDB()
553 {
554  // Set once per run
555  if ( fOADBSet ) return;
556 
557  AliDebug(1,"Begin");
558 
559  AliCDBManager *cdb = AliCDBManager::Instance();
560 
561  if (fOCDBpath.Length())
562  {
563  cdb->SetDefaultStorage(fOCDBpath.Data());
564  AliInfo(Form("Default storage %s",fOCDBpath.Data()));
565  }
566 
567  cdb->SetRun(fRun);
568 
569  //
570  // EMCAL from RAW OCDB
571  if (fOCDBpath.Contains("alien:"))
572  {
573  cdb->SetSpecificStorage("EMCAL/Calib/Data","raw://");
574  cdb->SetSpecificStorage("EMCAL/Calib/Time","raw://");
575  cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","raw://");
576  }
577 
578  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
579 
580  fOADBSet = kTRUE;
581 
582  return ;
583 }
584 
585 //_____________________________________________________
589 //_____________________________________________________
590 void AliAnalysisTaskEMCALClusterize::AddNewTCardInducedCellsToDigit()
591 {
592  // Get the number of stored digits, to assign the index of the new ones
593  Int_t idigit = fDigitsArr->GetEntriesFast();
594 
595  Double_t fixTime = 615.*1e-9;
596 
597  for(Int_t j = 0; j < fgkNEMCalCells; j++)
598  {
599  // Newly created?
600  if ( !fTCardCorrCellsNew[j] ) continue;
601 
602  // Accept only if at least 10 MeV
603  if ( fTCardCorrCellsEner[j] < 0.01 ) continue;
604 
605  // Check if it was not masked
606  Float_t amp = 0;
607  Double_t time = 0;
608  if ( !fRecoUtils->AcceptCalibrateCell(j,0,amp,time,fCaloCells) ) continue;
609 
610 // printf("add new digit absId %d, accept? %d, digit %d, induced amp %2.2f\n",
611 // j,fRecoUtils->AcceptCalibrateCell(j,0,amp,time,fCaloCells),idigit,fTCardCorrCellsEner[j]);
612 
613  // Now add the cell to the digits list
614  //AliEMCALDigit* digit =
615  new((*fDigitsArr)[idigit]) AliEMCALDigit( -2, -2, j, fTCardCorrCellsEner[j], fixTime,AliEMCALDigit::kHG,idigit, 0, 0, 0);
616 
617  //digit->SetListOfParents(0,0x0,0x0); // not needed
618 
619  idigit++;
620  }// loop on all possible cells
621 }
622 
623 
624 //_____________________________________________________
629 //_____________________________________________________
630 void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
631 {
632  fEvent = 0x0;
633 
634  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
635  Int_t eventN = Entry();
636  if(aodIH) eventN = aodIH->GetReadEntry();
637 
638  if (eventN > fMaxEvent)
639  return ;
640 
641  //printf("Clusterizer --- Event %d-- \n",eventN);
642 
643  // Check if input event are embedded events
644  // If so, take output event
645  if (aodIH && aodIH->GetMergeEvents())
646  {
647  fEvent = AODEvent();
648 
649  if(!aodIH->GetMergeEMCALCells())
650  AliFatal("Events merged but not EMCAL cells, check analysis settings!");
651 
652  AliDebug(1,"Use embedded events");
653 
654  AliDebug(1,Form("\t InputEvent N Clusters %d, N Cells %d",
655  InputEvent()->GetNumberOfCaloClusters(),InputEvent()->GetEMCALCells()->GetNumberOfCells()));
656 
657  AliDebug(1,Form("\t MergedEvent N Clusters %d, N Cells %d",
658  aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()));
659 
660 // if(DebugLevel() > 1)
661 // {
662 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
663 // {
664 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
665 // if(sigCluster->IsEMCAL()) AliInfo(Form("\t \t Signal cluster: i %d, E %f",icl,sigCluster->E()));
666 // }
667 // }
668 
669  AliDebug(1,Form("\t OutputEvent N Clusters %d, N Cells %d",
670  AODEvent()->GetNumberOfCaloClusters(), AODEvent()->GetEMCALCells()->GetNumberOfCells()));
671  }
672  else if(fInputFromFilter)
673  {
674  //printf("Get Input From Filtered AOD\n");
675  fEvent = AODEvent();
676  }
677  else
678  {
679  fEvent = InputEvent();
680  if(fFillAODCaloCells) FillAODCaloCells();
681  if(fFillAODHeader) FillAODHeader();
682  }
683 
684  if (!fEvent)
685  {
686  AliError("Event not available");
687  return ;
688  }
689 
690  //Recover the pointer to CaloCells container
691  if ( fInputCaloCellsName.Length() == 0 )
692  fCaloCells = fEvent->GetEMCALCells();
693  else
694  {
695  fCaloCells = (AliVCaloCells*) fEvent->FindListObject(fInputCaloCellsName);
696  if ( !fCaloCells )
697  AliWarning(Form("CaloCells branch <%s> not found use STD!",fInputCaloCellsName.Data()));
698  else
699  fCaloCells = fEvent->GetEMCALCells();
700  }
701 
702  //Process events if there is a high energy cluster
703  if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
704 
705  //-------------------------------------------------------------------------------------
706  // Reject events if LED was firing, use only for LHC11a data
707  // Reject event if triggered by exotic cell and remove exotic cells if not triggered
708  //-------------------------------------------------------------------------------------
709 
710  if( IsLEDEvent( fRun ) ) { fEvent = 0x0 ; return ; }
711 
712  if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
713 
714  //-------------------------------------------------------------------------------------
715  // Set the cluster array in the event (output or input)
716  //-------------------------------------------------------------------------------------
717 
718  if ( fFillAODFile )
719  {
720  //Magic line to write events to AOD filem put after event rejection
721  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
722  }
723  else if( !fOutputAODBranchSet )
724  {
725  // Create array of clusters/cells and put it in the input event, if output AOD not selected, only once
726  InputEvent()->AddObject(fOutputAODBranch);
727  AliInfo(Form("Add AOD clusters branch <%s> to input event",fOutputAODBranchName.Data()));
728 
729  if ( fOutputAODBranchName.Length() > 0 )
730  {
731  InputEvent()->AddObject(fOutputAODCells);
732  AliInfo(Form("Add AOD cells branch <%s> to input event",fOutputAODCellsName.Data()));
733  }
734 
735  fOutputAODBranchSet = kTRUE;
736  }
737 }
738 
739 //____________________________________________________
744 //____________________________________________________
745 void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
746 {
747  Int_t nClusters = fEvent->GetNumberOfCaloClusters();
748  Int_t nClustersOrg = 0;
749 
750  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
751  if(aodIH && aodIH->GetEventToMerge()) //Embedding
752  nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
753 
754  ResetArrays();
755 
756  // Loop on original clusters, get MC labels, cluster time (OLD AODs),
757  // or track matching residuals (if matching is not requested)
758  if(fSetCellMCLabelFromEdepFrac || fSetCellMCLabelFromCluster || fRecalibrateWithClusterTime || !fDoTrackMatching)
759  {
760  for (Int_t i = 0; i < nClusters; i++)
761  {
762  AliVCluster *clus = 0;
763  if(aodIH && aodIH->GetEventToMerge()) //Embedding
764  clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
765  else
766  clus = fEvent->GetCaloCluster(i);
767 
768  if(!clus) return;
769 
770  nClustersOrg++;
771 
772  if(!clus->IsEMCAL()) continue;
773 
774  Int_t label = clus->GetLabel();
775  Int_t label2 = -1 ;
776  if (clus->GetNLabels() >=2 ) label2 = clus->GetLabelAt(1) ;
777 
778  //printf("Org cluster %d) ID = %d, E %2.2f, Time %3.0f, N Cells %d, N MC labels %d, main MC label %d, all MC labels:\n",
779  // i, clus->GetID(), clus->E(), clus->GetTOF()*1e9, clus->GetNCells(), clus->GetNLabels(), label);
780  //
781  //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++)
782  // printf("%d) Label %d, E dep frac %0.2f; ",
783  // imc, clus->GetLabelAt(imc),clus->GetClusterMCEdepFraction(imc));
784  //if(clus->GetNLabels() > 0) printf("\n");
785 
786  UShort_t * index = clus->GetCellsAbsId() ;
787  for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
788  {
789  fOrgClusterCellId[index[icell]] = i;
790  fCellTime[index[icell]] = clus->GetTOF();
791  fCellMatchdEta[index[icell]] = clus->GetTrackDz();
792  fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
793 
794  if(!fSetCellMCLabelFromEdepFrac)
795  {
796  fCellLabels[index[icell]] = label;
797  fCellSecondLabels[index[icell]] = label2;
798  }
799  }
800  }
801  }
802 
803  // Do here induced cell energy assignation by T-Card correlation emulation, ONLY MC
804  if(fTCardCorrEmulation) MakeCellTCardCorrelation();
805 
806  //
807  // Transform CaloCells into Digits
808  //
809  Int_t idigit = 0;
810  Int_t id = -1;
811  Float_t amp = -1;
812  Double_t time = -1;
813 
814  Int_t bc = InputEvent()->GetBunchCrossNumber();
815 
816  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
817  {
818  // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
819  id = fCaloCells->GetCellNumber(icell);
820  Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,fCaloCells);
821  time-=fConstantTimeShift*1e-9; // only in case of simulations done before 2015
822 
823  // Do not include cells with too low energy, nor exotic cell
824  if( amp < fRecParam->GetMinECut() ||
825  time > fRecParam->GetTimeMax() ||
826  time < fRecParam->GetTimeMin() ) accept = kFALSE;
827 
828  // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
829  if (fRecalibrateWithClusterTime)
830  {
831  time = fCellTime[id];
832  //printf("cell %d time org %f - ",id, time*1.e9);
833  fRecoUtils->RecalibrateCellTime(id,bc,time);
834  //printf("recal %f\n",time*1.e9);
835  }
836 
837  //Exotic?
838  if (accept && fRecoUtils->IsExoticCell(id,fCaloCells,bc))
839  accept = kFALSE;
840 
841  if( !accept )
842  {
843  AliDebug(2,Form("Remove channel absId %d, index %d of %d, amp %f, time %f",
844  id,icell, fCaloCells->GetNumberOfCells(), amp, time*1.e9));
845  continue;
846  }
847 
848  //
849  // MC
850  //
851 
852  // Old way
853  Int_t mcLabel = fCaloCells->GetMCLabel(icell);
854  Float_t eDep = amp;
855 
856  //printf("--- Selected cell %d--- amp %f\n",id,amp);
857 
858  // New way
859  TArrayI labeArr(0);
860  TArrayF eDepArr(0);
861  Int_t nLabels = 0;
862 
863  if(!fSetCellMCLabelFromEdepFrac)
864  {
865  // Old way to recover/set the cell MC label
866  // Only possibility for old Run1 productions
867 
868  if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
869 
870  else if( fSetCellMCLabelFromCluster == 0 &&
871  fRemapMCLabelForAODs) RemapMCLabelForAODs(mcLabel);
872 
873  else mcLabel = -1; // found later
874 
875  // Last parameter of the digit object should be MC deposited energy,
876  // since it is not available in aliroot before year 2016, add just the cell amplitude so that
877  // we give more weight to the MC label of the cell with highest energy in the cluster
878 
879  Float_t efrac = fCaloCells->GetEFraction(icell);
880 
881  // When checking the MC of digits, give weight to cells with embedded signal
882  if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
883 
884  eDep *= efrac ;
885  }
886  else // fSetCellMCLabelFromEdepFrac = true
887  {
888  // New way, valid only for MC productions with aliroot > v5-07-21
889  mcLabel = -1;
890 
891  // Map the digit to cell index for later to calculate the cell MC energy deposition map
892  fCellLabels[id] = idigit;
893  //printf("\t absId %d, idigit %d\n",id,idigit);
894 
895  if(fOrgClusterCellId[id] < 0) continue; // index can be negative if noisy cell that did not form cluster
896 
897  AliVCluster *clus = 0;
898  Int_t iclus = fOrgClusterCellId[id];
899 
900  if(iclus < 0)
901  {
902  AliInfo("Negative original cluster index, skip \n");
903  continue;
904  }
905 
906  if(aodIH && aodIH->GetEventToMerge()) //Embedding
907  clus = aodIH->GetEventToMerge()->GetCaloCluster(iclus); //Get clusters directly from embedded signal
908  else
909  clus = fEvent->GetCaloCluster(iclus);
910 
911  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(id, clus, MCEvent(), amp, labeArr, eDepArr);
912  nLabels = labeArr.GetSize();
913 
914  } // cell MC label, new
915 
916  // Apply here the found induced energies
917  if ( fTCardCorrEmulation )
918  {
919 // if( TMath::Abs(fTCardCorrCellsEner[id]) > 0.001 )
920 // printf("add energy to digit %d, absId %d: amp %2.2f + %2.2f\n",idigit,id,amp,fTCardCorrCellsEner[id]);
921  amp+=fTCardCorrCellsEner[id];
922  }
923 
924  //
925  // Create the digit
926  //
927  if(amp <= 0.01) continue ; // accept if > 10 MeV
928 
929  //if(amp > 1.5)printf("*** Add digit *** amp %f, nlabels %d, label %d, found %d, edep fract tot %f, ncluslabels %d\n",amp, nLabels,mcLabel,found,edepTotFrac,nClusLabels);
930 
931  AliEMCALDigit* digit = new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, eDep);
932 
933  if(nLabels > 0)
934  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
935 
936  idigit++;
937 
938  }
939 
940  fDigitsArr->Sort();
941 
942  // Call after Sort, if not it screws up digits index order in clusterization
943  if ( fTCardCorrEmulation ) AddNewTCardInducedCellsToDigit();
944 
945 
946  //-------------------------------------------------------------------------------------
947  // Do the clusterization
948  //-------------------------------------------------------------------------------------
949 
950  fClusterizer->Digits2Clusters("");
951 
952  //-------------------------------------------------------------------------------------
953  // Transform the recpoints into AliVClusters
954  //-------------------------------------------------------------------------------------
955 
956  RecPoints2Clusters();
957 
958  if(!fCaloClusterArr)
959  {
960  AliWarning(Form("No array with CaloClusters, input RecPoints entries %d",fClusterArr->GetEntriesFast()));
961  return;
962  }
963 
964  AliDebug(1,Form("N clusters: before recluster %d, after recluster %d, recpoints %d",
965  nClustersOrg, fCaloClusterArr->GetEntriesFast(),fClusterArr->GetEntriesFast()));
966 
967 // if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
968 // {
969 // AliInfo("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
970 // fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
971 // fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
972 // }
973 }
974 
975 
976 //_____________________________________________________
978 //_____________________________________________________
979 void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
980 {
981  Double_t cellAmplitude = 0;
982  Double_t cellTime = 0;
983  Short_t cellNumber = 0;
984  Int_t cellMCLabel = 0;
985  Double_t cellEFrac = 0;
986  Int_t nClustersOrg = 0;
987 
988  // Fill the array with the EMCAL clusters, copy them
989  for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
990  {
991  AliVCluster *clus = fEvent->GetCaloCluster(i);
992  if(clus->IsEMCAL())
993  {
994  //recalibrate/remove bad channels/etc if requested
995  if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
996  {
997  continue;
998  }
999 
1000  if(fRecoUtils->IsRecalibrationOn())
1001  {
1002  //Calibrate cluster
1003  fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, fCaloCells);
1004 
1005  //CalibrateCells
1006  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1007  {
1008  if (fCaloCells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
1009  break;
1010 
1011  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1012  fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
1013  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1014 
1015  //Do not include bad channels found in analysis?
1016  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
1017  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
1018  continue;
1019 
1020  fCaloCells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
1021 
1022  }// cells loop
1023  }// recalibrate
1024 
1025  //Cast to ESD or AOD, needed to create the cluster array
1026  AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
1027  AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
1028 
1029  if (esdCluster)
1030  {
1031  fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
1032  }//ESD
1033  else if(aodCluster)
1034  {
1035  fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
1036  }//AOD
1037  else
1038  AliWarning("Wrong CaloCluster type?");
1039 
1040  nClustersOrg++;
1041  }
1042  }
1043 
1044  //Do the unfolding
1045  fUnfolder->UnfoldClusters(fCaloClusterArr, fCaloCells);
1046 
1047  //CLEAN-UP
1048  fUnfolder->Clear();
1049 }
1050 
1051 //_______________________________________________________________
1064 //_______________________________________________________________
1065 void AliAnalysisTaskEMCALClusterize::ConfigureEMCALRecoUtils
1066 (Bool_t bMC , Bool_t bExotic, Bool_t bNonLin,
1067  Bool_t bRecalE, Bool_t bBad , Bool_t bRecalT, Int_t debug)
1068 {
1069  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - **** Start ***\n");
1070 
1071  // Init
1072  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils ;
1073 
1074  // Exotic cells removal
1075 
1076  if(bExotic)
1077  {
1078  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Remove exotics in EMCAL\n");
1079  fRecoUtils->SwitchOnRejectExoticCell() ;
1080  fRecoUtils->SwitchOnRejectExoticCluster();
1081 
1082  // fRecoUtils->SetExoticCellDiffTimeCut(50); // If |t cell max - t cell in cross| > 50 do not add its energy, avoid
1083  fRecoUtils->SetExoticCellFractionCut(0.97); // 1-Ecross/Ecell > 0.97 -> out
1084  fRecoUtils->SetExoticCellMinAmplitudeCut(4.); // 4 GeV
1085  }
1086 
1087  // Recalibration factors
1088 
1089  if(bRecalE && ! bMC)
1090  {
1091  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on energy recalibration in EMCAL\n");
1092  fRecoUtils->SwitchOnRecalibration();
1093  fRecoUtils->SwitchOnRunDepCorrection();
1094  }
1095 
1096  // Remove EMCAL hot channels
1097 
1098  if(bBad)
1099  {
1100  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on bad channels removal in EMCAL\n");
1101  fRecoUtils->SwitchOnBadChannelsRemoval();
1102  fRecoUtils->SwitchOnDistToBadChannelRecalculation();
1103  }
1104 
1105  // *** Time recalibration settings ***
1106 
1107  if(bRecalT && ! bMC)
1108  {
1109  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() - Switch on time recalibration in EMCAL\n");
1110  fRecoUtils->SwitchOnTimeRecalibration();
1111  fRecoUtils->SwitchOnL1PhaseInTimeRecalibration() ;
1112  }
1113 
1114  // Recalculate position with method
1115 
1116  fRecoUtils->SetPositionAlgorithm(AliEMCALRecoUtils::kPosTowerGlobal);
1117 
1118  // Non linearity
1119 
1120  if( bNonLin )
1121  {
1122  if(!bMC)
1123  {
1124  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kBeamTestCorrected xxx\n");
1125  fRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrectedv3);
1126  }
1127  else
1128  {
1129  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx SET Non linearity correction kPi0MCv3 xxx\n");
1130  fRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kPi0MCv3);
1131  }
1132  }
1133  else
1134  {
1135  if ( debug > 0 ) printf("ConfigureEMCALRecoUtils() xxx DON'T SET Non linearity correction xxx\n");
1136  fRecoUtils->SetNonLinearityFunction(AliEMCALRecoUtils::kNoCorrection);
1137  }
1138 
1139 }
1140 
1141 //_____________________________________________________
1143 //_____________________________________________________
1144 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells()
1145 {
1146  AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
1147  Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
1148 
1149  AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
1150  aodEMcells.CreateContainer(nEMcell);
1151  aodEMcells.SetType(AliVCaloCells::kEMCALCell);
1152  Double_t calibFactor = 1.;
1153  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1154  {
1155  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1156  fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
1157  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1158 
1159  if(fRecoUtils->IsRecalibrationOn())
1160  {
1161  calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
1162  }
1163 
1164  if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
1165  { //Channel is not declared as bad
1166  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
1167  eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
1168  }
1169  else
1170  {
1171  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
1172  }
1173  }
1174 
1175  aodEMcells.Sort();
1176 }
1177 
1178 //__________________________________________________
1180 //__________________________________________________
1181 void AliAnalysisTaskEMCALClusterize::FillAODHeader()
1182 {
1183  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
1184  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
1185 
1186  Double_t pos[3] ;
1187  Double_t covVtx[6];
1188  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
1189 
1190  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
1191  if(!header) AliFatal("Not a standard AOD");
1192  header->SetRunNumber(fRun);
1193 
1194  if(esdevent)
1195  {
1196  TTree* tree = fInputHandler->GetTree();
1197  if (tree)
1198  {
1199  TFile* file = tree->GetCurrentFile();
1200  if (file) header->SetESDFileName(file->GetName());
1201  }
1202  }
1203  else if (aodevent) {
1204  AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
1205  if(!aodheader) AliFatal("Not a standard AOD");
1206  header->SetESDFileName(aodheader->GetESDFileName());
1207  }
1208 
1209  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
1210  header->SetOrbitNumber(fEvent->GetOrbitNumber());
1211  header->SetPeriodNumber(fEvent->GetPeriodNumber());
1212  header->SetEventType(fEvent->GetEventType());
1213 
1214  // Centrality
1215  if(fUseAliCentrality)
1216  {
1217  if(fEvent->GetCentrality())
1218  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
1219  else
1220  header->SetCentrality(0);
1221  }
1222 
1223  // Trigger
1224  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
1225 
1226  header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
1227 
1228  header->SetTriggerMask(fEvent->GetTriggerMask());
1229  header->SetTriggerCluster(fEvent->GetTriggerCluster());
1230 
1231  if (esdevent)
1232  {
1233  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
1234  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
1235  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
1236  }
1237  else if (aodevent)
1238  {
1239  header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
1240  header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
1241  header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
1242  }
1243 
1244  header->SetMagneticField(fEvent->GetMagneticField());
1245  //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
1246 
1247  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
1248  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
1249  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
1250  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
1251  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
1252 
1253  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
1254  Float_t diamcov[3];
1255  fEvent->GetDiamondCovXY(diamcov);
1256  header->SetDiamond(diamxy,diamcov);
1257  if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
1258  else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
1259 
1260  //
1261  Int_t nVertices = 1 ;/* = prim. vtx*/;
1262  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1263 
1264  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1265 
1266  // Access to the AOD container of vertices
1267  TClonesArray &vertices = *(AODEvent()->GetVertices());
1268  Int_t jVertices=0;
1269 
1270  // Add primary vertex. The primary tracks will be defined
1271  // after the loops on the composite objects (V0, cascades, kinks)
1272  fEvent->GetPrimaryVertex()->GetXYZ(pos);
1273  Float_t chi = 0;
1274  if (esdevent)
1275  {
1276  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1277  chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
1278  }
1279  else if (aodevent)
1280  {
1281  aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1282  chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
1283  }
1284 
1285  AliAODVertex * primary = new(vertices[jVertices++])
1286  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
1287  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
1288  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
1289 }
1290 
1291 //___________________________________________________________
1295 //___________________________________________________________
1296 void AliAnalysisTaskEMCALClusterize::FillCaloClusterInEvent()
1297 {
1298  // First recalculate track-matching for the new clusters
1299  if(fDoTrackMatching)
1300  {
1301  fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom,MCEvent());
1302  }
1303 
1304  // Put the new clusters in the AOD list
1305 
1306  Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
1307 
1308  for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
1309  {
1310  AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1311 
1312  newCluster->SetID(i);
1313 
1314  // Correct cluster energy non linearity
1315  newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
1316 
1317  //Add matched track
1318  if(fDoTrackMatching)
1319  {
1320  Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1321  if(trackIndex >= 0)
1322  {
1323  newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1324  AliDebug(2,Form("Matched Track index %d to new cluster %d",trackIndex,i));
1325  }
1326 
1327  Float_t dR = 999., dZ = 999.;
1328  fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dZ,dR);
1329  newCluster->SetTrackDistance(dR,dZ);
1330  }
1331  else
1332  {// Assign previously assigned matched track in reco, very very rough
1333  Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1334  newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1335  }
1336 
1337  //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1338  //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1339  //printf("\n");
1340  //
1341  //printf("New cluster %d) ID = %d, E %2.2f, Time %3.0f, N Cells %d, N MC labels %d, main MC label %d, all MC labels:\n",
1342  // i, newCluster->GetID(), newCluster->E(), newCluster->GetTOF()*1e9, newCluster->GetNCells(), newCluster->GetNLabels(), newCluster->GetLabel());
1343  //
1344  //for(Int_t imc = 0; imc < newCluster->GetNLabels(); imc++)
1345  // printf("%d) Label %d, E dep frac %0.2f; ",
1346  // imc,newCluster->GetLabelAt(imc),newCluster->GetClusterMCEdepFraction(imc));
1347  //if(newCluster->GetNLabels() > 0) printf("\n");
1348 
1349  // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1350  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, newCluster);
1351 
1352  new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
1353 
1354  AliDebug(2,Form("New cluster %d of %d, energy %f, mc label %d",
1355  newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()));
1356 
1357  } // cluster loop
1358 
1359  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1360 }
1361 
1362 
1363 //________________________________________________________________
1365 //________________________________________________________________
1366 Float_t AliAnalysisTaskEMCALClusterize::GetEventCentrality() const
1367 {
1368  if(fUseAliCentrality)
1369  {
1370  if(GetCentrality()) return GetCentrality()->GetCentralityPercentile(fCentralityClass) ;
1371  else return -1. ;
1372  }
1373  else
1374  {
1375  if(GetMultSelCen()) return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass,kTRUE) ;
1376  else return -1. ;
1377  }
1378 
1379  return -1.;
1380 }
1381 
1382 //_______________________________________________
1384 //_______________________________________________
1385 TString AliAnalysisTaskEMCALClusterize::GetPass()
1386 {
1387  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1388  {
1389  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null");
1390  return TString("");
1391  }
1392 
1393  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1394  {
1395  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null");
1396  return TString("");
1397  }
1398 
1399  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1400  if (pass.Contains("ass1")) return TString("pass1");
1401  else if (pass.Contains("ass2")) return TString("pass2");
1402  else if (pass.Contains("ass3")) return TString("pass3");
1403  else if (pass.Contains("ass4")) return TString("pass4");
1404  else if (pass.Contains("ass5")) return TString("pass5");
1405  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1406  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1407  {
1408  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1409  return TString("pass1");
1410  }
1411  else if (pass.Contains("LHC14a1a"))
1412  {
1413  AliInfo("Enable EMCal energy calibration for this MC production!!");
1414 
1415  fRecoUtils->SwitchOnRecalibration();
1416 
1417  return TString("LHC14a1a");
1418  }
1419 
1420  // No condition fullfilled, give a default value
1421  AliInfo("Pass number string not found");
1422 
1423  return TString("");
1424 }
1425 
1426 //_________________________________________
1429 //_________________________________________
1430 void AliAnalysisTaskEMCALClusterize::Init()
1431 {
1432  if(fDebug >=0) (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1433 
1434  fOADBSet = kFALSE;
1435  if(fOADBFilePath == "") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
1436 
1437  fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1438 
1439  if(!fRecParam) fRecParam = new AliEMCALRecParam;
1440  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1441 
1442  if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1443  if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1444  if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1445  fRejectBelowThreshold = kFALSE;
1446 
1447  // Centrality
1448  if(fCentralityClass == "") fCentralityClass = "V0M";
1449 
1450  if (fOCDBpath == "") fOCDBpath = "raw://" ;
1451  if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1452 
1453  if(gROOT->LoadMacro(fConfigName) >=0)
1454  {
1455  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
1456  AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1457  fGeomName = clus->fGeomName;
1458  fLoadGeomMatrices = clus->fLoadGeomMatrices;
1459  fOCDBpath = clus->fOCDBpath;
1460  fAccessOCDB = clus->fAccessOCDB;
1461  fRecParam = clus->fRecParam;
1462  fJustUnfold = clus->fJustUnfold;
1463  fFillAODFile = clus->fFillAODFile;
1464  fRecoUtils = clus->fRecoUtils;
1465  fConfigName = clus->fConfigName;
1466  fMaxEvent = clus->fMaxEvent;
1467  fDoTrackMatching = clus->fDoTrackMatching;
1468  fUpdateCell = clus->fUpdateCell;
1469  fOutputAODBranchName = clus->fOutputAODBranchName;
1470  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1471  fCentralityClass = clus->fCentralityClass;
1472  fCentralityBin[0] = clus->fCentralityBin[0];
1473  fCentralityBin[1] = clus->fCentralityBin[1];
1474  fUseAliCentrality = clus->fUseAliCentrality;
1475  }
1476 }
1477 
1478 //_______________________________________________________
1481 //_______________________________________________________
1482 void AliAnalysisTaskEMCALClusterize::InitClusterization()
1483 {
1484  if (fJustUnfold)
1485  {
1486  // init the unfolding afterburner
1487  delete fUnfolder;
1488  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1489  return;
1490  }
1491 
1492  //First init the clusterizer
1493  delete fClusterizer;
1494  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1495  fClusterizer = new AliEMCALClusterizerv1 (fGeom);
1496  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1497  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1498  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1499  {
1500  fClusterizer = new AliEMCALClusterizerNxN(fGeom);
1501  fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1502  fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1503  }
1504  else
1505  {
1506  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1507  }
1508 
1509  //Now set the parameters
1510  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1511  fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1512  fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1513  fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1514  fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1515  fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1516  fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1517  fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1518  fClusterizer->SetTimeCalibration ( fRecParam->IsTimeCalibrationOn() );
1519  fClusterizer->SetInputCalibrated ( kTRUE );
1520  fClusterizer->SetJustClusters ( kTRUE );
1521 
1522  // Initialize the cluster rec points and digits arrays and get them.
1523  fClusterizer->SetOutput(0);
1524  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1525  fDigitsArr = fClusterizer->GetDigits();
1526 
1527  //In case of unfolding after clusterization is requested, set the corresponding parameters
1528  if(fRecParam->GetUnfold())
1529  {
1530  Int_t i=0;
1531  for (i = 0; i < 8; i++)
1532  {
1533  fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1534  }//end of loop over parameters
1535 
1536  for (i = 0; i < 3; i++)
1537  {
1538  fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1539  fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1540  }//end of loop over parameters
1541 
1542  fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1543  fClusterizer->InitClusterUnfolding();
1544 
1545  }// to unfold
1546 }
1547 
1548 //________________________________________________________________
1552 //________________________________________________________________
1553 void AliAnalysisTaskEMCALClusterize::InitGeometry()
1554 {
1555  if(fGeomMatrixSet) return;
1556 
1557  if (!fGeom)
1558  {
1559  if(fGeomName=="")
1560  {
1561  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(fRun);
1562  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fGeom->GetName(),fRun));
1563  }
1564  else
1565  {
1566  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1567  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeomName.Data()));
1568  }
1569 
1570  // Init geometry, I do not like much to do it like this ...
1571  if(fImportGeometryFromFile && !gGeoManager)
1572  {
1573  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1574  {
1575  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1576  if (fRun < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1577  else if (fRun < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1578  else if (fRun < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1579  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >=2015
1580  }
1581 
1582  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1583 
1584  TGeoManager::Import(fImportGeometryFilePath) ;
1585  }
1586 
1587  AliDebug(1,Form("Init for run=%d",fRun));
1588  if (!gGeoManager) AliDebug(1,"Careful!, gGeoManager not loaded, load misalign matrices");
1589  } // geometry pointer did not exist before
1590 
1591  if(fLoadGeomMatrices)
1592  {
1593  AliInfo("Load user defined EMCAL geometry matrices");
1594 
1595  // OADB if available
1596  AliOADBContainer emcGeoMat("AliEMCALgeo");
1597  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1598  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(fRun,"EmcalMatrices");
1599 
1600  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1601  {
1602 
1603  if (!fGeomMatrix[mod]) // Get it from OADB
1604  {
1605  AliDebug(2,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
1606  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1607 
1608  fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1609  }
1610 
1611  if(fGeomMatrix[mod])
1612  {
1613  if(DebugLevel() > 1)
1614  fGeomMatrix[mod]->Print();
1615 
1616  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1617  }
1618  else if(gGeoManager)
1619  {
1620  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
1621  fGeom->SetMisalMatrix(fGeom->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
1622  }
1623  else
1624  {
1625  AliError(Form("Alignment atrix for SM %d is not available",mod));
1626  }
1627 
1628  fGeomMatrixSet=kTRUE;
1629 
1630  }//SM loop
1631  }//Load matrices
1632  else if(!gGeoManager)
1633  {
1634  AliInfo("Get geo matrices from data");
1635  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1636  if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1637  {
1638  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
1639  }//AOD
1640  else
1641  {
1642  AliDebug(1,"Load Misaligned matrices");
1643 
1644  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1645 
1646  if(!esd)
1647  {
1648  AliError("This event does not contain ESDs?");
1649  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1650  return;
1651  }
1652 
1653  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1654  {
1655  if(DebugLevel() > 1)
1656  esd->GetEMCALMatrix(mod)->Print();
1657 
1658  if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1659 
1660  }
1661 
1662  fGeomMatrixSet=kTRUE;
1663 
1664  } // ESD
1665  } // Load matrices from Data
1666 }
1667 
1668 //____________________________________________________
1673 //____________________________________________________
1674 Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1675 {
1676  if(!fRemoveExoticEvents) return kFALSE;
1677 
1678  // Loop on cells
1679 
1680  Float_t totCellE = 0;
1681  Int_t bc = InputEvent()->GetBunchCrossNumber();
1682 
1683  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1684  {
1685  Float_t ecell = 0 ;
1686  Double_t tcell = 0 ;
1687 
1688  Int_t absID = fCaloCells->GetCellNumber(icell);
1689  Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,fCaloCells);
1690  tcell-=fConstantTimeShift*1e-9;// Only for MC simulations done before 2015
1691 
1692  if(accept && !fRecoUtils->IsExoticCell(absID,fCaloCells,bc)) totCellE += ecell;
1693  }
1694 
1695  // TString triggerclasses = event->GetFiredTriggerClasses();
1696  // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1697  // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1698  // return;
1699  //
1700 
1701  //printf("TotE cell %f\n",totCellE);
1702  if(totCellE < 1) return kTRUE;
1703 
1704  return kFALSE;
1705 }
1706 
1707 //________________________________________________________________
1709 //________________________________________________________________
1710 Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
1711 {
1712  if(!fRemoveLEDEvents) return kFALSE;
1713 
1714  //check events of LHC11a period
1715  if(run < 146858 || run > 146860) return kFALSE ;
1716 
1717  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1718  Int_t ncellsSM3 = 0;
1719  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1720  {
1721  if ( fCaloCells->GetAmplitude (icell) > 0.1 &&
1722  fCaloCells->GetCellNumber(icell)/(24*48)==3 ) ncellsSM3++;
1723  }
1724 
1725  TString triggerclasses = fEvent->GetFiredTriggerClasses();
1726 
1727  Int_t ncellcut = 21;
1728  if(triggerclasses.Contains("EMC")) ncellcut = 35;
1729 
1730  if( ncellsSM3 >= ncellcut)
1731  {
1732  AliInfo(Form("Reject event %d with ncells in SM3 %d",(Int_t)Entry(),ncellsSM3));
1733  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1734  return kTRUE;
1735  }
1736 
1737  return kFALSE;
1738 }
1739 
1740 
1741 //_______________________________________________________
1744 //_______________________________________________________
1745 void AliAnalysisTaskEMCALClusterize::MakeCellTCardCorrelation()
1746 {
1747  Int_t id = -1;
1748  Float_t amp = -1;
1749 
1750  // Loop on all cells with signal
1751  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1752  {
1753  id = fCaloCells->GetCellNumber(icell);
1754  amp = fCaloCells->GetAmplitude (icell); // fCaloCells->GetCellAmplitude(id);
1755 
1756  if ( amp <= fTCardCorrMinAmp ) continue ;
1757 
1758  //
1759  // First get the SM, col-row of this tower
1760  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1761  fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
1762  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1763 
1764  //
1765  // Determine randomly if we want to create a correlation for this cell,
1766  // depending the SM number of the cell
1767  Float_t rand = fRandom.Uniform(0, 1);
1768 
1769  if ( rand > fTCardCorrInduceEnerProb[imod] ) continue;
1770 
1771  AliDebug(1,Form("Reference cell absId %d, iEta %d, iPhi %d, amp %2.3f",id,ieta,iphi,amp));
1772 
1773  //
1774  // Get the absId of the cells in the cross and same T-Card
1775  Int_t absIDup = -1;
1776  Int_t absIDdo = -1;
1777  Int_t absIDlr = -1;
1778  Int_t absIDuplr = -1;
1779  Int_t absIDdolr = -1;
1780 
1781  Int_t absIDup2 = -1;
1782  Int_t absIDup2lr = -1;
1783  Int_t absIDdo2 = -1;
1784  Int_t absIDdo2lr = -1;
1785 
1786  // Only 2 columns in the T-Card, +1 for even and -1 for odd with respect reference cell
1787  Int_t colShift = 0;
1788  if ( (ieta%2) && ieta <= AliEMCALGeoParams::fgkEMCALCols-1 ) colShift = -1;
1789  if ( !(ieta%2) && ieta >= 0 ) colShift = +1;
1790 
1791  absIDlr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+colShift);
1792 
1793  // Check up / down cells from reference cell not out of SM and in same T-Card
1794  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1 )
1795  {
1796  absIDup = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
1797  absIDuplr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta+colShift);
1798  }
1799 
1800  if ( iphi > 0 )
1801  {
1802  absIDdo = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
1803  absIDdolr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta+colShift);
1804  }
1805 
1806  // Check 2 up / 2 down cells from reference cell not out of SM
1807  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-2 )
1808  {
1809  absIDup2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta);
1810  absIDup2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta+colShift);
1811  }
1812 
1813  if ( iphi > 1 )
1814  {
1815  absIDdo2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta);
1816  absIDdo2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta+colShift);
1817  }
1818 
1819  // In same T-Card?
1820  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+1)/8) ) { absIDup = -1 ; absIDuplr = -1 ; }
1821  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-1)/8) ) { absIDdo = -1 ; absIDdolr = -1 ; }
1822  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+2)/8) ) { absIDup2 = -1 ; absIDup2lr = -1 ; }
1823  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-2)/8) ) { absIDdo2 = -1 ; absIDdo2lr = -1 ; }
1824 
1825  //
1826  // Check if they are not declared bad or exist
1827  Bool_t okup = AcceptCell(absIDup );
1828  Bool_t okdo = AcceptCell(absIDdo );
1829  Bool_t oklr = AcceptCell(absIDlr );
1830  Bool_t okuplr = AcceptCell(absIDuplr );
1831  Bool_t okdolr = AcceptCell(absIDdolr );
1832  Bool_t okup2 = AcceptCell(absIDup2 );
1833  Bool_t okdo2 = AcceptCell(absIDdo2 );
1834  Bool_t okup2lr= AcceptCell(absIDup2lr);
1835  Bool_t okdo2lr= AcceptCell(absIDdo2lr);
1836 
1837  AliDebug(1,Form("Same T-Card cells:\n \t up %d (%d), down %d (%d), left-right %d (%d), up-lr %d (%d), down-lr %d (%d)\n"
1838  "\t up2 %d (%d), down2 %d (%d), up2-lr %d (%d), down2-lr %d (%d)",
1839  absIDup ,okup ,absIDdo ,okdo ,absIDlr,oklr,absIDuplr ,okuplr ,absIDdolr ,okdolr ,
1840  absIDup2,okup2,absIDdo2,okdo2, absIDup2lr,okup2lr,absIDdo2lr,okdo2lr));
1841 
1842  //
1843  // Generate some energy for the nearby cells in same TCard , depending on this cell energy
1844  // Check if originally the tower had no or little energy, in which case tag it as new
1845  Float_t fracupdown = fTCardCorrInduceEnerFrac[0][imod]+amp*fTCardCorrInduceEnerFracP1[0][imod];
1846  Float_t fracupdownleri = fTCardCorrInduceEnerFrac[1][imod]+amp*fTCardCorrInduceEnerFracP1[1][imod];
1847  Float_t fracleri = fTCardCorrInduceEnerFrac[2][imod]+amp*fTCardCorrInduceEnerFracP1[2][imod];
1848  Float_t frac2nd = fTCardCorrInduceEnerFrac[3][imod]+amp*fTCardCorrInduceEnerFracP1[3][imod];
1849 
1850  AliDebug(1,Form("Fraction for SM %d (min %2.3f, max %2.3f):\n"
1851  "\t up-down : c %2.3e, p1 %2.3e, p2 %2.4e, sig %2.3e, fraction %2.3f\n"
1852  "\t up-down-lr: c %2.3e, p1 %2.3e, p2 %2.4e, sig %2.3e, fraction %2.3f\n"
1853  "\t left-right: c %2.3e, p1 %2.3e, p2 %2.4e, sig %2.3e, fraction %2.3f\n"
1854  "\t 2nd row : c %2.3e, p1 %2.3e, p2 %2.4e, sig %2.3e, fraction %2.3f",
1855  imod, fTCardCorrInduceEnerFracMin[imod], fTCardCorrInduceEnerFracMax[imod],
1856  fTCardCorrInduceEner[0][imod],fTCardCorrInduceEnerFrac[0][imod],fTCardCorrInduceEnerFracP1[0][imod],fTCardCorrInduceEnerFracWidth[0][imod],fracupdown,
1857  fTCardCorrInduceEner[1][imod],fTCardCorrInduceEnerFrac[1][imod],fTCardCorrInduceEnerFracP1[1][imod],fTCardCorrInduceEnerFracWidth[1][imod],fracupdownleri,
1858  fTCardCorrInduceEner[2][imod],fTCardCorrInduceEnerFrac[2][imod],fTCardCorrInduceEnerFracP1[2][imod],fTCardCorrInduceEnerFracWidth[2][imod],fracleri,
1859  fTCardCorrInduceEner[3][imod],fTCardCorrInduceEnerFrac[3][imod],fTCardCorrInduceEnerFracP1[3][imod],fTCardCorrInduceEnerFracWidth[3][imod],frac2nd));
1860 
1861  if( fracupdown < fTCardCorrInduceEnerFracMin[imod] ) fracupdown = fTCardCorrInduceEnerFracMin[imod];
1862  if( fracupdown > fTCardCorrInduceEnerFracMax[imod] ) fracupdown = fTCardCorrInduceEnerFracMax[imod];
1863  if( fracupdownleri < fTCardCorrInduceEnerFracMin[imod] ) fracupdownleri = fTCardCorrInduceEnerFracMin[imod];
1864  if( fracupdownleri > fTCardCorrInduceEnerFracMax[imod] ) fracupdownleri = fTCardCorrInduceEnerFracMax[imod];
1865  if( fracleri < fTCardCorrInduceEnerFracMin[imod] ) fracleri = fTCardCorrInduceEnerFracMin[imod];
1866  if( fracleri > fTCardCorrInduceEnerFracMax[imod] ) fracleri = fTCardCorrInduceEnerFracMax[imod];
1867  if( frac2nd < fTCardCorrInduceEnerFracMin[imod] ) frac2nd = fTCardCorrInduceEnerFracMin[imod];
1868  if( frac2nd > fTCardCorrInduceEnerFracMax[imod] ) frac2nd = fTCardCorrInduceEnerFracMax[imod];
1869 
1870  // Randomize the induced fraction, if requested
1871  if(fRandomizeTCard)
1872  {
1873  fracupdown = fRandom.Gaus(fracupdown ,fTCardCorrInduceEnerFracWidth[0][imod]);
1874  fracupdownleri = fRandom.Gaus(fracupdownleri,fTCardCorrInduceEnerFracWidth[1][imod]);
1875  fracleri = fRandom.Gaus(fracleri ,fTCardCorrInduceEnerFracWidth[2][imod]);
1876  frac2nd = fRandom.Gaus(frac2nd ,fTCardCorrInduceEnerFracWidth[3][imod]);
1877 
1878  AliDebug(1,Form("Randomized fraction: up-down %2.3f; up-down-left-right %2.3f; left-right %2.3f; 2nd row %2.3f",
1879  fracupdown,fracupdownleri,fracleri,frac2nd));
1880  }
1881 
1882  // Calculate induced energy
1883  Float_t indEupdown = fTCardCorrInduceEner[0][imod]+amp*fracupdown;
1884  Float_t indEupdownleri = fTCardCorrInduceEner[1][imod]+amp*fracupdownleri;
1885  Float_t indEleri = fTCardCorrInduceEner[2][imod]+amp*fracleri;
1886  Float_t indE2nd = fTCardCorrInduceEner[3][imod]+amp*frac2nd;
1887 
1888  AliDebug(1,Form("Induced energy: up-down %2.3f; up-down-left-right %2.3f; left-right %2.3f; 2nd row %2.3f",
1889  indEupdown,indEupdownleri,indEleri,indE2nd));
1890 
1891  // Check if we induce too much energy, in such case use a constant value
1892  if ( fTCardCorrMaxInduced < indE2nd ) indE2nd = fTCardCorrMaxInduced;
1893  if ( fTCardCorrMaxInduced < indEupdownleri ) indEupdownleri = fTCardCorrMaxInduced;
1894  if ( fTCardCorrMaxInduced < indEupdown ) indEupdown = fTCardCorrMaxInduced;
1895  if ( fTCardCorrMaxInduced < indEleri ) indEleri = fTCardCorrMaxInduced;
1896 
1897  AliDebug(1,Form("Induced energy, saturated?: up-down %2.3f; up-down-left-right %2.3f; left-right %2.3f; 2nd row %2.3f",
1898  indEupdown,indEupdownleri,indEleri,indE2nd));
1899 
1900  //
1901  // Add the induced energy, check if cell existed
1902  if ( okup )
1903  {
1904  fTCardCorrCellsEner[absIDup] += indEupdown;
1905 
1906  if ( fCaloCells->GetCellAmplitude(absIDup) < 0.01 ) fTCardCorrCellsNew[absIDup] = kTRUE;
1907  }
1908 
1909  if ( okdo )
1910  {
1911  fTCardCorrCellsEner[absIDdo] += indEupdown;
1912 
1913  if ( fCaloCells->GetCellAmplitude(absIDdo) < 0.01 ) fTCardCorrCellsNew[absIDdo] = kTRUE;
1914  }
1915 
1916  if ( oklr )
1917  {
1918  fTCardCorrCellsEner[absIDlr] += indEleri;
1919 
1920  if ( fCaloCells->GetCellAmplitude(absIDlr) < 0.01 ) fTCardCorrCellsNew[absIDlr] = kTRUE;
1921  }
1922 
1923  if ( okuplr )
1924  {
1925  fTCardCorrCellsEner[absIDuplr] += indEupdownleri;
1926 
1927  if ( fCaloCells->GetCellAmplitude(absIDuplr ) < 0.01 ) fTCardCorrCellsNew[absIDuplr] = kTRUE;
1928  }
1929 
1930  if ( okdolr )
1931  {
1932  fTCardCorrCellsEner[absIDdolr] += indEupdownleri;
1933 
1934  if ( fCaloCells->GetCellAmplitude(absIDdolr ) < 0.01 ) fTCardCorrCellsNew[absIDdolr] = kTRUE;
1935  }
1936 
1937  if ( okup2 )
1938  {
1939  fTCardCorrCellsEner[absIDup2] += indE2nd;
1940 
1941  if ( fCaloCells->GetCellAmplitude(absIDup2) < 0.01 ) fTCardCorrCellsNew[absIDup2] = kTRUE;
1942  }
1943 
1944  if ( okup2lr )
1945  {
1946  fTCardCorrCellsEner[absIDup2lr] += indE2nd;
1947 
1948  if ( fCaloCells->GetCellAmplitude(absIDup2lr) < 0.01 ) fTCardCorrCellsNew[absIDup2lr] = kTRUE;
1949  }
1950 
1951  if ( okdo2 )
1952  {
1953  fTCardCorrCellsEner[absIDdo2] += indE2nd;
1954 
1955  if ( fCaloCells->GetCellAmplitude(absIDdo2) < 0.01 ) fTCardCorrCellsNew[absIDdo2] = kTRUE;
1956  }
1957 
1958  if ( okdo2lr )
1959  {
1960  fTCardCorrCellsEner[absIDdo2lr] += indE2nd;
1961 
1962  if ( fCaloCells->GetCellAmplitude(absIDdo2lr) < 0.01 ) fTCardCorrCellsNew[absIDdo2lr] = kTRUE;
1963  }
1964 
1965  //
1966  // Subtract the added energy to main cell, if energy conservation is requested
1967  if ( fTCardCorrClusEnerConserv )
1968  {
1969  if ( oklr ) fTCardCorrCellsEner[id] -= indEleri;
1970  if ( okuplr ) fTCardCorrCellsEner[id] -= indEupdownleri;
1971  if ( okdolr ) fTCardCorrCellsEner[id] -= indEupdownleri;
1972  if ( okup ) fTCardCorrCellsEner[id] -= indEupdown;
1973  if ( okdo ) fTCardCorrCellsEner[id] -= indEupdown;
1974  if ( okup2 ) fTCardCorrCellsEner[id] -= indE2nd;
1975  if ( okup2lr ) fTCardCorrCellsEner[id] -= indE2nd;
1976  if ( okdo2 ) fTCardCorrCellsEner[id] -= indE2nd;
1977  if ( okdo2lr ) fTCardCorrCellsEner[id] -= indE2nd;
1978  } // conserve energy
1979 
1980  } // cell loop
1981 
1982 }
1983 
1984 //_______________________________________________________
1986 //_______________________________________________________
1987 void AliAnalysisTaskEMCALClusterize::PrintParam()
1988 {
1989  AliInfo(Form("Geometry: name <%s>, matrix set <%d>, load matrix <%d>, import geo <%d> from path <%s>",
1990  fGeomName.Data(), fGeomMatrixSet, fLoadGeomMatrices, fImportGeometryFromFile, fImportGeometryFilePath.Data()));
1991 
1992  if ( fAccessOCDB ) AliInfo(Form("OCDB path name <%s>", fOCDBpath.Data()));
1993  if ( fAccessOADB ) AliInfo(Form("OADB path name <%s>", fOADBFilePath.Data()));
1994 
1995  if ( fInputCaloCellsName.Length() > 0 )
1996  AliInfo(Form("Input CaloCells <%s>", fInputCaloCellsName.Data()));
1997 
1998  AliInfo(Form("Just Unfold clusters <%d>, new clusters list name <%s>, new cells name <%s>",
1999  fJustUnfold, fOutputAODBranchName.Data(), fOutputAODCellsName.Data()));
2000 
2001  if ( fFillAODFile ) AliInfo(Form("Fill new AOD file with: header <%d>, cells <%d>",fFillAODHeader,fFillAODCaloCells));
2002 
2003  AliInfo(Form("Use cell time for cluster <%d>, Apply constant time shift <%2.2f>, Do track-matching <%d>, Update cells <%d>, Input from ESD filter <%d>",
2004  fRecalibrateWithClusterTime, fConstantTimeShift, fDoTrackMatching, fUpdateCell, fInputFromFilter));
2005 
2006  AliInfo(Form("Reject events: larger than <%d>, LED <%d>, exotics <%d>", fMaxEvent, fRemoveLEDEvents, fRemoveExoticEvents));
2007 
2008  if (fCentralityBin[0] != -1 && fCentralityBin[1] != -1 )
2009  AliInfo(Form("Centrality bin [%2.2f,%2.2f], class <%s>, use AliCentrality? <%d>",
2010  fCentralityBin[0], fCentralityBin[1], fCentralityClass.Data(), fUseAliCentrality));
2011 
2012  if ( fSelectEMCALEvent )
2013  AliInfo(Form("Select events with signal in EMCal: E min <%2.2f>, n cell min <%d>", fEMCALEnergyCut, fEMCALNcellsCut));
2014 
2015  AliInfo(Form("MC label from cluster <%d>, Use EdepFrac <%d>, remap AODs <%d>",
2016  fSetCellMCLabelFromCluster, fSetCellMCLabelFromEdepFrac, fRemapMCLabelForAODs));
2017 }
2018 
2019 //_______________________________________________________
2021 //_______________________________________________________
2022 void AliAnalysisTaskEMCALClusterize::PrintTCardParam()
2023 {
2024  if(!fTCardCorrEmulation)
2025  {
2026  AliInfo("T-Card emulation not activated");
2027  return;
2028  }
2029 
2030  AliInfo(Form("T-Card emulation activated, energy conservation <%d>, randomize E <%d>, induced energy parameters:",
2031  fTCardCorrClusEnerConserv,fRandomizeTCard));
2032 
2033  AliInfo(Form("T-Card emulation super-modules fraction: Min cell E %2.2f Max induced E %2.2f",
2034  fTCardCorrMinAmp,fTCardCorrMaxInduced));
2035 
2036  for(Int_t ism = 0; ism < 22; ism++)
2037  {
2038  printf("\t sm %d, fraction %2.3f, E frac abs min %2.3e max %2.3e \n",
2039  ism, fTCardCorrInduceEnerProb[ism],fTCardCorrInduceEnerFracMin[ism],fTCardCorrInduceEnerFracMax[ism]);
2040 
2041  for(Int_t icell = 0; icell < 4; icell++)
2042  {
2043  printf("\t \t cell type %d, c %2.4e, p0 %2.4e, p1 %2.4e, sigma %2.4e \n",
2044  icell,fTCardCorrInduceEner[icell][ism],fTCardCorrInduceEnerFrac[icell][ism],
2045  fTCardCorrInduceEnerFracP1[icell][ism],fTCardCorrInduceEnerFracWidth[icell][ism]);
2046  }
2047  }
2048 }
2049 
2050 //_______________________________________________________
2054 //_______________________________________________________
2055 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters()
2056 {
2057  Int_t j = 0;
2058  for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
2059  {
2060  AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
2061 
2062  const Int_t ncells = recPoint->GetMultiplicity();
2063  Int_t ncellsTrue = 0;
2064 
2065  if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
2066 
2067  // cells and their amplitude fractions
2068  UShort_t absIds[ncells];
2069  Double32_t ratios[ncells];
2070 
2071  //For later check embedding
2072  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2073 
2074  Float_t clusterE = 0;
2075  for (Int_t c = 0; c < ncells; c++)
2076  {
2077  AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
2078 
2079  absIds[ncellsTrue] = digit->GetId();
2080  ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
2081 
2082  if ( !fRecParam->GetUnfold() && (ratios[ncellsTrue] > 1 || ratios[ncellsTrue] < 1) )
2083  AliWarning(Form("recpoint cell E %2.3f but digit E %2.3f and no unfolding", recPoint->GetEnergiesList()[c], digit->GetAmplitude()));
2084 
2085  // In case of unfolding, remove digits with unfolded energy too low
2086  if(fSelectCell)
2087  {
2088  if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
2089  {
2090 
2091  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
2092  recPoint->GetEnergiesList()[c],digit->GetAmplitude()));
2093  continue;
2094  } // if cuts
2095  }// Select cells
2096 
2097  //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
2098  clusterE +=recPoint->GetEnergiesList()[c];
2099 
2100  // In case of embedding, fill ratio with amount of signal,
2101  if (aodIH && aodIH->GetMergeEvents())
2102  {
2103  //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
2104  AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
2105  AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
2106 
2107  Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2108  //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2109  Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2110  //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
2111 
2112  if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
2113  //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
2114 
2115  }//Embedding
2116 
2117  ncellsTrue++;
2118 
2119  }// cluster cell loop
2120 
2121  if (ncellsTrue < 1)
2122  {
2123  AliDebug(2,Form("Skipping cluster with no cells avobe threshold E = %f, ncells %d",
2124  recPoint->GetEnergy(), ncells));
2125  continue;
2126  }
2127 
2128  //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
2129 
2130  if(clusterE < fRecParam->GetClusteringThreshold())
2131  {
2132  AliDebug(2,Form("Remove cluster with energy below seed threshold %f",clusterE));
2133  continue;
2134  }
2135 
2136  TVector3 gpos;
2137  Float_t g[3];
2138 
2139  // calculate new cluster position
2140 
2141  recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
2142  recPoint->GetGlobalPosition(gpos);
2143  gpos.GetXYZ(g);
2144 
2145  // create a new cluster
2146 
2147  (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
2148  AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
2149  j++;
2150  clus->SetType(AliVCluster::kEMCALClusterv1);
2151  clus->SetE(clusterE);
2152  clus->SetPosition(g);
2153  clus->SetNCells(ncellsTrue);
2154  clus->SetCellsAbsId(absIds);
2155  clus->SetCellsAmplitudeFraction(ratios);
2156  clus->SetChi2(-1); //not yet implemented
2157  clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
2158  clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
2159  clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
2160 
2161  if(ncells == ncellsTrue)
2162  {
2163  Float_t elipAxis[2];
2164  recPoint->GetElipsAxis(elipAxis);
2165  clus->SetM02(elipAxis[0]*elipAxis[0]) ;
2166  clus->SetM20(elipAxis[1]*elipAxis[1]) ;
2167  clus->SetDispersion(recPoint->GetDispersion());
2168  }
2169  else if(fSelectCell)
2170  {
2171  // In case some cells rejected, in unfolding case, recalculate
2172  // shower shape parameters and position
2173  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
2174 
2175  AliVCaloCells* cells = 0x0;
2176  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
2177  else cells = InputEvent()->GetEMCALCells();
2178 
2179  fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
2180  fRecoUtils->RecalculateClusterPID(clus);
2181  fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
2182  }
2183 
2184  // MC
2185 
2186  if ( fSetCellMCLabelFromCluster == 1 ) SetClustersMCLabelFrom2SelectedLabels(recPoint,clus) ;
2187  else if( fSetCellMCLabelFromCluster == 2 ) SetClustersMCLabelFromOriginalClusters(clus) ;
2188  else
2189  {
2190  //
2191  // Normal case, trust what the clusterizer has found
2192  //
2193  Int_t parentMult = 0;
2194  Int_t *parentList = recPoint->GetParents(parentMult);
2195  Float_t *parentListDE = recPoint->GetParentsDE(); // deposited energy
2196 
2197  clus->SetLabel(parentList, parentMult);
2198  if(fSetCellMCLabelFromEdepFrac)
2199  clus->SetClusterMCEdepFractionFromEdepArray(parentListDE);
2200 
2201  //
2202  // Set the cell energy deposition fraction map:
2203  //
2204  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
2205  {
2206  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
2207 
2208  // Get the digit that originated this cell cluster
2209 // AliVCaloCells* cells = 0x0;
2210 // if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
2211 // else cells = InputEvent()->GetEMCALCells();
2212 
2213  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
2214  {
2215  Int_t idigit = fCellLabels[absIds[icell]];
2216 
2217  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
2218 
2219  // Find the 4 MC labels that contributed to the cluster and their
2220  // deposited energy in the current digit
2221 
2222  mcEdepFracPerCell[icell] = 0; // init
2223 
2224  Int_t nparents = dig->GetNiparent();
2225  if ( nparents > 0 )
2226  {
2227  Int_t digLabel =-1 ;
2228  Float_t edep = 0 ;
2229  Float_t edepTot = 0 ;
2230  Float_t mcEDepFrac[4] = {0,0,0,0};
2231 
2232  // all parents in digit
2233  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2234  {
2235  digLabel = dig->GetIparent (jndex+1);
2236  edep = dig->GetDEParent(jndex+1);
2237  edepTot += edep;
2238 
2239  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
2240  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
2241  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
2242  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
2243  } // all prarents in digit
2244 
2245  // Divide energy deposit by total deposited energy
2246  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
2247  if(edepTot > 0.01)
2248  {
2249  mcEdepFracPerCell[icell] = clus->PackMCEdepFraction(mcEDepFrac);
2250  }
2251  } // at least one parent label in digit
2252  } // cell in cluster loop
2253 
2254  clus->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
2255 
2256  delete [] mcEdepFracPerCell;
2257 
2258  } // at least one parent in cluster, do the cell primary packing
2259  }
2260  } // recPoints loop
2261 }
2262 
2263 //_____________________________________________________________________
2266 //_____________________________________________________________________
2267 void AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs(Int_t & label)
2268 {
2269  if(label < 0) return ;
2270 
2271  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
2272  if(!evt) return ;
2273 
2274  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
2275  if(!arr) return ;
2276 
2277  if(label < arr->GetEntriesFast())
2278  {
2279  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
2280  if(!particle) return ;
2281 
2282  if(label == particle->Label()) return ; // label already OK
2283  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
2284  }
2285  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
2286 
2287  // loop on the particles list and check if there is one with the same label
2288  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
2289  {
2290  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
2291  if(!particle) continue ;
2292 
2293  if(label == particle->Label())
2294  {
2295  label = ind;
2296  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
2297  return;
2298  }
2299  }
2300 
2301  label = -1;
2302 
2303  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
2304 }
2305 
2306 //________________________________________________
2308 //________________________________________________
2309 void AliAnalysisTaskEMCALClusterize::ResetArrays()
2310 {
2311  for(Int_t j = 0; j < fgkNEMCalCells; j++)
2312  {
2313  fOrgClusterCellId[j] = -1;
2314  fCellLabels[j] = -1;
2315  fCellSecondLabels[j] = -1;
2316  fCellTime[j] = 0.;
2317  fCellMatchdEta[j] = -999;
2318  fCellMatchdPhi[j] = -999;
2319 
2320  fTCardCorrCellsEner[j] = 0.;
2321  fTCardCorrCellsNew [j] = kFALSE;
2322  }
2323 }
2324 
2325 //_____________________________________________________________________________________________________
2329 //_____________________________________________________________________________________________________
2330 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint * recPoint,
2331  AliAODCaloCluster * clus)
2332 {
2333  Int_t parentMult = 0;
2334  Int_t *parentList = recPoint->GetParents(parentMult);
2335  clus->SetLabel(parentList, parentMult);
2336 
2337  //Write the second major contributor to each MC cluster.
2338  Int_t iNewLabel ;
2339  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2340  {
2341 
2342  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2343  if(idCell>=0)
2344  {
2345  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
2346  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
2347  {
2348  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
2349  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
2350  }
2351  if (iNewLabel == 1)
2352  {
2353  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
2354  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
2355  {
2356  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
2357  }
2358 
2359  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
2360  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
2361  delete [] newLabelArray;
2362  }
2363  } // positive cell number
2364  }
2365 }
2366 
2367 //___________________________________________________________________________________________________
2371 //___________________________________________________________________________________________________
2372 void AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster * clus)
2373 {
2374  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
2375  clArray.Reset();
2376  Int_t nClu = 0;
2377  Int_t nLabTotOrg = 0;
2378  Float_t emax = -1;
2379  Int_t idMax = -1;
2380 
2381  AliVEvent * event = fEvent;
2382 
2383  // In case of embedding the MC signal is in the added event
2384  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2385  if(aodIH && aodIH->GetEventToMerge()) //Embedding
2386  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
2387 
2388 
2389  // Find the clusters that originally had the cells
2390  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2391  {
2392  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2393 
2394  if(idCell>=0)
2395  {
2396  Int_t idCluster = fOrgClusterCellId[idCell];
2397 
2398  Bool_t set = kTRUE;
2399  for(Int_t icl =0; icl < nClu; icl++)
2400  {
2401  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
2402  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
2403  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
2404  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
2405  }
2406  if( set && idCluster >= 0)
2407  {
2408  clArray.SetAt(idCluster,nClu++);
2409  //printf("******** idCluster %d \n",idCluster);
2410  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
2411 
2412  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
2413 
2414  //Search highest E cluster
2415  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2416  //printf("\t E %f\n",clOrg->E());
2417  if(emax < clOrg->E())
2418  {
2419  emax = clOrg->E();
2420  idMax = idCluster;
2421  }
2422  }
2423  }
2424  }// cell loop
2425 
2426  if(nClu==0 || nLabTotOrg == 0)
2427  {
2428  //if(clus->E() > 0.25) printf("AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters() - Check: N org clusters %d, n tot labels %d, cluster E %f, n cells %d\n",nClu,nLabTotOrg,clus->E(), clus->GetNCells());
2429  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
2430  //printf("\n");
2431  }
2432 
2433  // Put the first in the list the cluster with highest energy
2434  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
2435  {
2436  Int_t maxIndex = -1;
2437  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
2438  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2439  {
2440  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
2441  }
2442 
2443  if(firstCluster >=0 && idMax >=0)
2444  {
2445  clArray.SetAt(idMax,0);
2446  clArray.SetAt(firstCluster,maxIndex);
2447  }
2448  }
2449 
2450  // Get the labels list in the original clusters, assign all to the new cluster
2451  TArrayI clMCArray(nLabTotOrg) ;
2452  clMCArray.Reset();
2453 
2454  Int_t nLabTot = 0;
2455  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2456  {
2457  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
2458  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
2459  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2460  Int_t nLab = clOrg->GetNLabels();
2461 
2462  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
2463  {
2464  Int_t lab = clOrg->GetLabelAt(iLab) ;
2465  if(lab>=0)
2466  {
2467  Bool_t set = kTRUE;
2468  //printf("\t \t Set Label %d \n", lab);
2469  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
2470  {
2471  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
2472  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
2473  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
2474  }
2475  if( set ) clMCArray.SetAt(lab,nLabTot++);
2476  }
2477  }
2478  }// cluster loop
2479 
2480  // Set the final list of labels
2481 
2482  clus->SetLabel(clMCArray.GetArray(), nLabTot);
2483 
2484 // printf("Final list of labels for new cluster : \n");
2485 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
2486 // {
2487 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
2488 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
2489 // printf(" org %d ",label);
2490 // RemapMCLabelForAODs(label);
2491 // printf(" new %d \n",label);
2492 // }
2493 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
2494 }
2495 
2496 //____________________________________________________________
2497 // Init geometry, create list of output clusters and cells.
2498 //____________________________________________________________
2499 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
2500 {
2501  // Clusters
2502  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2503 
2504  if(fOutputAODBranchName.Length()==0)
2505  {
2506  fOutputAODBranchName = "newEMCALClustersArray";
2507  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
2508  }
2509 
2510  fOutputAODBranch->SetName(fOutputAODBranchName);
2511 
2512  if( fFillAODFile )
2513  {
2514  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2515 
2516  //fOutputAODBranch->SetOwner(kFALSE);
2517 
2518  AddAODBranch("TClonesArray", &fOutputAODBranch);
2519  }
2520 
2521  // Cells
2522  if ( fOutputAODBranchName.Length() > 0 )
2523  {
2524  fOutputAODCells = new AliAODCaloCells(fOutputAODCellsName,fOutputAODCellsName,AliAODCaloCells::kEMCALCell);
2525 
2526  if( fFillAODFile ) AddAODBranch("AliAODCaloCells", &fOutputAODCells);
2527  }
2528 }
2529 
2530 //_______________________________________________________________________
2533 //________________________________________________________________________
2534 void AliAnalysisTaskEMCALClusterize::UpdateCells()
2535 {
2536  // Update cells only in case re-calibration was done
2537  // or bad map applied or additional T-Card cells added.
2538  if(!fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
2539  !fRecoUtils->IsRecalibrationOn() &&
2540  !fRecoUtils->IsRunDepRecalibrationOn() &&
2541  !fRecoUtils->IsTimeRecalibrationOn() &&
2542  !fRecoUtils->IsL1PhaseInTimeRecalibrationOn() &&
2543  !fTCardCorrEmulation ) return;
2544 
2545  const Int_t ncells = fCaloCells->GetNumberOfCells();
2546  const Int_t ndigis = fDigitsArr->GetEntries();
2547 
2548  if ( fOutputAODCellsName.Length() > 0 )
2549  {
2550  fOutputAODCells->DeleteContainer();
2551  fOutputAODCells->CreateContainer(ndigis);
2552  }
2553  else if ( ncells != ndigis ) // update case
2554  {
2555  fCaloCells->DeleteContainer();
2556  fCaloCells->CreateContainer(ndigis);
2557  }
2558 
2559  for (Int_t idigit = 0; idigit < ndigis; ++idigit)
2560  {
2561  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
2562 
2563  Double_t cellAmplitude = digit->GetAmplitude();
2564  Short_t cellNumber = digit->GetId();
2565  Double_t cellTime = digit->GetTime();
2566 
2567  Bool_t highGain = kFALSE;
2568  if( digit->GetType() == AliEMCALDigit::kHG ) highGain = kTRUE;
2569 
2570  // Only for MC
2571  // Get the label of the primary particle that generated the cell
2572  // Assign the particle that deposited more energy
2573  Int_t nparents = digit->GetNiparent();
2574  Int_t cellMcEDepFrac =-1 ;
2575  Float_t cellMcLabel =-1.;
2576  if ( nparents > 0 )
2577  {
2578  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2579  {
2580  if(cellMcEDepFrac >= digit->GetDEParent(jndex+1)) continue ;
2581 
2582  cellMcLabel = digit->GetIparent (jndex+1);
2583  cellMcEDepFrac= digit->GetDEParent(jndex+1);
2584  } // all primaries in digit
2585  } // select primary label
2586 
2587  if ( cellMcEDepFrac < 0 ) cellMcEDepFrac = 0.;
2588 
2589  if ( fUpdateCell )
2590  fCaloCells ->SetCell(idigit, cellNumber, cellAmplitude, cellTime, cellMcLabel, cellMcEDepFrac, highGain);
2591  else
2592  fOutputAODCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime, cellMcLabel, cellMcEDepFrac, highGain);
2593  }
2594 }
2595 
2596 //_______________________________________________________
2607 //_______________________________________________________
2608 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
2609 {
2610  //-------
2611  // Step 1
2612 
2613  // Remove the contents of AOD branch output list set in the previous event
2614  fOutputAODBranch->Clear("C");
2615 
2616  LoadBranches();
2617 
2618  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
2619  // If we need a centrality bin, we select only those events in the corresponding bin.
2620  if( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
2621  {
2622  Float_t cen = GetEventCentrality();
2623  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
2624  }
2625 
2626  // Intermediate array with new clusters : init the array only once or clear from previous event
2627  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
2628  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
2629 
2630 
2631  // In case of analysis in multiple runs, check the OADB again
2632  if ( InputEvent()->GetRunNumber() != fRun )
2633  {
2634  fRun = InputEvent()->GetRunNumber();
2635 
2636  fOADBSet = kFALSE; // recover the OADB for this run
2637 
2638  AliInfo(Form("Set run to %d",fRun));
2639  }
2640 
2641  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
2642 
2643  // Get the event, do some checks and settings
2644  CheckAndGetEvent() ;
2645 
2646  if (!fEvent)
2647  {
2648  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
2649  return ;
2650  }
2651 
2652  // Init pointers, geometry, clusterizer, ocdb, aodb
2653 
2654  if(fAccessOCDB) AccessOCDB();
2655  if(fAccessOADB) AccessOADB(); // only once
2656 
2657  InitClusterization();
2658 
2659  // Print once the analysis parameters
2660  if ( fDebug > 0 || !fPrintOnce )
2661  {
2662  fRecParam->Print("reco");
2663 
2664  PrintParam();
2665 
2666  PrintTCardParam();
2667 
2668  fPrintOnce = kTRUE;
2669  }
2670 
2671  //-------
2672  // Step 2
2673 
2674  // Make clusters
2675  if (fJustUnfold) ClusterUnfolding();
2676  else ClusterizeCells() ;
2677 
2678  //-------
2679  // Step 3
2680 
2681  FillCaloClusterInEvent();
2682 
2683  if ( fUpdateCell || fOutputAODCellsName.Length() > 0 )
2684  UpdateCells();
2685 
2686 }
2687 
2688 
AliAnalysisTaskEMCALClusterize::fOutputAODCellsName
TString fOutputAODCellsName
New of output cells AOD branch name.
Definition: AliAnalysisTaskEMCALClusterize.h:353
AliAnalysisTaskEMCALClusterize::FillAODCaloCells
void FillAODCaloCells()
Put calo cells in standard branch.
Definition: AliAnalysisTaskEMCALClusterize.cxx:1144
AliAnalysisTaskEMCALClusterize::fTCardCorrMinAmp
Float_t fTCardCorrMinAmp
Minimum cell energy to induce signal on adjacent cells.
Definition: AliAnalysisTaskEMCALClusterize.h:439
AliAnalysisTaskEMCALClusterize::fCaloCells
AliVCaloCells * fCaloCells
! CaloCells container
Definition: AliAnalysisTaskEMCALClusterize.h:341
AliAnalysisTaskEMCALClusterize::fTCardCorrEmulation
Bool_t fTCardCorrEmulation
Activate T-Card cells energy correlation.
Definition: AliAnalysisTaskEMCALClusterize.h:422
AliAnalysisTaskEMCALClusterize::fTCardCorrInduceEnerFracP1
Float_t fTCardCorrInduceEnerFracP1[4][22]
Induced energy loss gauss fraction param1 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
Definition: AliAnalysisTaskEMCALClusterize.h:429
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskEMCALClusterize::ResetArrays
virtual void ResetArrays()
Reset arrays containing information for all possible cells.
Definition: AliAnalysisTaskEMCALClusterize.cxx:2309
TH2F
Definition: External.C:236
debug
const int debug
Definition: scanAll.C:15
AliAnalysisTaskEMCALClusterize::GetPass
TString GetPass()
Get or guess pass number/string from path of filename.
Definition: AliAnalysisTaskEMCALClusterize.cxx:1385
AliAnalysisTaskEMCALClusterize.h
AliAnalysisTaskEMCALClusterize::fFillAODHeader
Bool_t fFillAODHeader
Copy header to standard branch.
Definition: AliAnalysisTaskEMCALClusterize.h:358
AliAnalysisTaskEMCALClusterize::fDigitsArr
TClonesArray * fDigitsArr
! Digits array
Definition: AliAnalysisTaskEMCALClusterize.h:338
AliAnalysisTaskEMCALClusterize::fTCardCorrClusEnerConserv
Bool_t fTCardCorrClusEnerConserv
When making correlation, subtract from the reference cell the induced energy on the neighbour cells...
Definition: AliAnalysisTaskEMCALClusterize.h:423
AliAnalysisTaskEMCALClusterize::fClusterizer
AliEMCALClusterizer * fClusterizer
! EMCAL clusterizer
Definition: AliAnalysisTaskEMCALClusterize.h:345
AliAnalysisTaskEMCALClusterize::fJustUnfold
Bool_t fJustUnfold
Just unfold, do not recluster.
Definition: AliAnalysisTaskEMCALClusterize.h:347
AliAnalysisTaskEMCALClusterize::fImportGeometryFromFile
Bool_t fImportGeometryFromFile
Import geometry settings in geometry.root file.
Definition: AliAnalysisTaskEMCALClusterize.h:388
AliAnalysisTaskEMCALClusterize::fSelectCell
Bool_t fSelectCell
Reject cells from cluster if energy is too low and recalculate position/energy and other...
Definition: AliAnalysisTaskEMCALClusterize.h:381
AliAnalysisTaskEMCALClusterize::fConfigName
TString fConfigName
Name of analysis configuration file.
Definition: AliAnalysisTaskEMCALClusterize.h:364
AliAnalysisTaskEMCALClusterize::fTCardCorrInduceEnerFracMax
Float_t fTCardCorrInduceEnerFracMax[22]
In case fTCardCorrInduceEnerFracP1 is non null, restrict the maximum fraction of induced energy per S...
Definition: AliAnalysisTaskEMCALClusterize.h:431
c
TCanvas * c
Definition: TestFitELoss.C:172
AliAnalysisTaskEMCALClusterize::PrintTCardParam
void PrintTCardParam()
Print parameters for T-Card correlation emulation.
Definition: AliAnalysisTaskEMCALClusterize.cxx:2022
AliAnalysisTaskEMCALClusterize::fCellSecondLabels
Int_t fCellSecondLabels[fgkNEMCalCells]
Array with Second MC label to be passed to digit.
Definition: AliAnalysisTaskEMCALClusterize.h:370
AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFrom2SelectedLabels
void SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint *recPoint, AliAODCaloCluster *clus)
Definition: AliAnalysisTaskEMCALClusterize.cxx:2330
AliAnalysisTaskEMCALClusterize::fCellMatchdEta
Float_t fCellMatchdEta[fgkNEMCalCells]
Array with cluster-track dPhi.
Definition: AliAnalysisTaskEMCALClusterize.h:372
AliAnalysisTaskEMCALClusterize::fOutputAODBranch
TClonesArray * fOutputAODBranch
! AOD Branch with output clusters
Definition: AliAnalysisTaskEMCALClusterize.h:350
AliAnalysisTaskEMCALClusterize::fOADBFilePath
TString fOADBFilePath
Default path $ALICE_PHYSICS/OADB/EMCAL, if needed change.
Definition: AliAnalysisTaskEMCALClusterize.h:393
AliAnalysisTaskEMCALClusterize::fTCardCorrCellsEner
Float_t fTCardCorrCellsEner[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells.
Definition: AliAnalysisTaskEMCALClusterize.h:424
AliAnalysisTaskEMCALClusterize::fImportGeometryFilePath
TString fImportGeometryFilePath
path fo geometry.root file
Definition: AliAnalysisTaskEMCALClusterize.h:389
AliAnalysisTaskEMCALClusterize::fTCardCorrInduceEnerFrac
Float_t fTCardCorrInduceEnerFrac[4][22]
Induced energy loss gauss fraction param0 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
Definition: AliAnalysisTaskEMCALClusterize.h:428
AliAnalysisTaskEMCALClusterize::ConfigureEMCALRecoUtils
void ConfigureEMCALRecoUtils(Bool_t bMC=kFALSE, Bool_t bExotic=kTRUE, Bool_t bNonLin=kFALSE, Bool_t bRecalE=kTRUE, Bool_t bBad=kTRUE, Bool_t bRecalT=kTRUE, Int_t debug=-1)
Definition: AliAnalysisTaskEMCALClusterize.cxx:1066
AliAnalysisTaskEMCALClusterize::fPrintOnce
Bool_t fPrintOnce
Print once analysis parameters.
Definition: AliAnalysisTaskEMCALClusterize.h:442
AliAnalysisTaskEMCALClusterize::fSelectEMCALEvent
Bool_t fSelectEMCALEvent
Process the event if there is some high energy cluster.
Definition: AliAnalysisTaskEMCALClusterize.h:402
AliAnalysisTaskEMCALClusterize::fTCardCorrMaxInduced
Float_t fTCardCorrMaxInduced
Maximum induced energy signal on adjacent cells.
Definition: AliAnalysisTaskEMCALClusterize.h:440
AliAnalysisTaskEMCALClusterize::GetEventCentrality
Float_t GetEventCentrality() const
Get centrality/multiplicity percentile.
Definition: AliAnalysisTaskEMCALClusterize.cxx:1366
AliAnalysisTaskEMCALClusterize::fInputFromFilter
Bool_t fInputFromFilter
Get the input from AODs from the filter.
Definition: AliAnalysisTaskEMCALClusterize.h:418
AliAnalysisTaskEMCALClusterize::fOutputAODBranchName
TString fOutputAODBranchName
New of output clusters AOD branch.
Definition: AliAnalysisTaskEMCALClusterize.h:351
AliAnalysisTaskEMCALClusterize::fGeom
AliEMCALGeometry * fGeom
EMCAL geometry.
Definition: AliAnalysisTaskEMCALClusterize.h:327
AliAnalysisTaskEMCALClusterize::ClusterUnfolding
void ClusterUnfolding()
Take the event clusters and unfold them.
Definition: AliAnalysisTaskEMCALClusterize.cxx:979
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskEMCALClusterize::fAccessOADB
Bool_t fAccessOADB
Get calibration from OADB for EMCAL.
Definition: AliAnalysisTaskEMCALClusterize.h:392
AliAnalysisTaskEMCALClusterize::fRemoveLEDEvents
Bool_t fRemoveLEDEvents
Remove LED events, use only for LHC11a.
Definition: AliAnalysisTaskEMCALClusterize.h:385
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliAnalysisTaskEMCALClusterize::fRecoUtils
AliEMCALRecoUtils * fRecoUtils
Access to factorized reconstruction algorithms.
Definition: AliAnalysisTaskEMCALClusterize.h:363
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskEMCALClusterize::fRejectBelowThreshold
Bool_t fRejectBelowThreshold
split (false-default) or reject (true) cell energy below threshold after UF
Definition: AliAnalysisTaskEMCALClusterize.h:384
AliAnalysisTaskEMCALClusterize::fConstantTimeShift
Float_t fConstantTimeShift
Apply a 600 ns time shift in case of simulation, shift in ns.
Definition: AliAnalysisTaskEMCALClusterize.h:394
AliAnalysisTaskEMCALClusterize::SetClustersMCLabelFromOriginalClusters
void SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster *clus)
Definition: AliAnalysisTaskEMCALClusterize.cxx:2372
AliAnalysisTaskEMCALClusterize::fUpdateCell
Bool_t fUpdateCell
On/Off the upate of the CaloCells container.
Definition: AliAnalysisTaskEMCALClusterize.h:380
AliAnalysisTaskEMCALClusterize::fSelectCellMinE
Float_t fSelectCellMinE
Min energy cell threshold, after unfolding.
Definition: AliAnalysisTaskEMCALClusterize.h:382
AliAnalysisTaskEMCALClusterize::fOADBSet
Bool_t fOADBSet
AODB parameters already set.
Definition: AliAnalysisTaskEMCALClusterize.h:391
AliAnalysisTaskEMCALClusterize::fOutputAODCells
AliAODCaloCells * fOutputAODCells
! AOD Branch with output cells
Definition: AliAnalysisTaskEMCALClusterize.h:352
AliAnalysisTaskEMCALClusterize::fTCardCorrInduceEner
Float_t fTCardCorrInduceEner[4][22]
Induced energy loss gauss constant on 0-same row, diff col, 1-up/down cells left/right col 2-left/rig...
Definition: AliAnalysisTaskEMCALClusterize.h:427
AliAnalysisTaskEMCALClusterize::fCellLabels
Int_t fCellLabels[fgkNEMCalCells]
Array with MC label to be passed to digit.
Definition: AliAnalysisTaskEMCALClusterize.h:369
AliAnalysisTaskEMCALClusterize::fRemoveExoticEvents
Bool_t fRemoveExoticEvents
Remove exotic events.
Definition: AliAnalysisTaskEMCALClusterize.h:386
AliAnalysisTaskEMCALClusterize::fDoTrackMatching
Bool_t fDoTrackMatching
On/Off the matching recalculation to speed up analysis in PbPb.
Definition: AliAnalysisTaskEMCALClusterize.h:379
AliAnalysisTaskEMCALClusterize::fCellMatchdPhi
Float_t fCellMatchdPhi[fgkNEMCalCells]
Array with cluster-track dEta.
Definition: AliAnalysisTaskEMCALClusterize.h:373
AliAnalysisTaskEMCALClusterize::fTCardCorrInduceEnerProb
Float_t fTCardCorrInduceEnerProb[22]
Probability to induce energy loss per SM.
Definition: AliAnalysisTaskEMCALClusterize.h:433
AliAnalysisTaskEMCALClusterize::fInputCaloCellsName
TString fInputCaloCellsName
Input cells branch name, if different from default branch.
Definition: AliAnalysisTaskEMCALClusterize.h:354
AliAnalysisTaskEMCALClusterize::fGeomMatrix
TGeoHMatrix * fGeomMatrix[22]
Geometry matrices with alignments.
Definition: AliAnalysisTaskEMCALClusterize.h:329
AliAnalysisTaskEMCALClusterize::IsLEDEvent
Bool_t IsLEDEvent(const Int_t run)
Check if event is LED, is so remove it. Affected LHC11a runs.
Definition: AliAnalysisTaskEMCALClusterize.cxx:1710
AliAnalysisTaskEMCALClusterize::fRecParam
AliEMCALRecParam * fRecParam
Reconstruction parameters container.
Definition: AliAnalysisTaskEMCALClusterize.h:344
AliAnalysisTaskEMCALClusterize::fGeomName
TString fGeomName
Name of geometry to use.
Definition: AliAnalysisTaskEMCALClusterize.h:328
AliAnalysisTaskEMCALClusterize::fOutputAODBranchSet
Bool_t fOutputAODBranchSet
Set the AOD clusters branch in the input event once.
Definition: AliAnalysisTaskEMCALClusterize.h:355
AliAnalysisTaskEMCALClusterize::PrintParam
void PrintParam()
Print clusterization task parameters.
Definition: AliAnalysisTaskEMCALClusterize.cxx:1987
Short_t
short Short_t
Definition: External.C:23
AliAnalysisTaskEMCALClusterize::fOCDBpath
TString fOCDBpath
Path with OCDB location.
Definition: AliAnalysisTaskEMCALClusterize.h:334
AliAnalysisTaskEMCALClusterize::fUnfolder
AliEMCALAfterBurnerUF * fUnfolder
! Unfolding procedure
Definition: AliAnalysisTaskEMCALClusterize.h:346
AliAnalysisTaskEMCALClusterize::fGeomMatrixSet
Bool_t fGeomMatrixSet
Set geometry matrices only once, for the first event.
Definition: AliAnalysisTaskEMCALClusterize.h:330
AliAnalysisTaskEMCALClusterize::fCentralityBin
Float_t fCentralityBin[2]
Minimum and maximum value of the centrality for the analysis.
Definition: AliAnalysisTaskEMCALClusterize.h:398
AliAnalysisTaskEMCALClusterize::fFillAODCaloCells
Bool_t fFillAODCaloCells
Copy calocells to standard branch.
Definition: AliAnalysisTaskEMCALClusterize.h:359
AliAnalysisTaskEMCALClusterize::fCaloClusterArr
TObjArray * fCaloClusterArr
! CaloClusters array
Definition: AliAnalysisTaskEMCALClusterize.h:340
AliAnalysisTaskEMCALClusterize::fCellTime
Double_t fCellTime[fgkNEMCalCells]
Array with cluster time to be passed to digit in case of AODs.
Definition: AliAnalysisTaskEMCALClusterize.h:371
AliAnalysisTaskEMCALClusterize::fCentralityClass
TString fCentralityClass
Name of selected centrality class.
Definition: AliAnalysisTaskEMCALClusterize.h:397
file
TFile * file
TList with histograms for a given trigger.
Definition: DrawAnaCaloTrackQA.C:89
AliAnalysisTaskEMCALClusterize::fAccessOCDB
Bool_t fAccessOCDB
Need to access info from OCDB (not really)
Definition: AliAnalysisTaskEMCALClusterize.h:335
AliAnalysisTaskEMCALClusterize::fClusterArr
TObjArray * fClusterArr
! Recpoints array
Definition: AliAnalysisTaskEMCALClusterize.h:339
UShort_t
unsigned short UShort_t
Definition: External.C:28
AliAnalysisTaskEMCALClusterize::fRecalibrateWithClusterTime
Bool_t fRecalibrateWithClusterTime
Use fCellTime to store time of cells in cluster.
Definition: AliAnalysisTaskEMCALClusterize.h:375
GetRunNumber
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
Option_t
const char Option_t
Definition: External.C:48
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskEMCALClusterize::FillAODHeader
void FillAODHeader()
Put event header information in standard AOD branch.
Definition: AliAnalysisTaskEMCALClusterize.cxx:1181
AliAnalysisTaskEMCALClusterize::fgkNEMCalCells
static const Int_t fgkNEMCalCells
Total number of cells in the calorimeter, 10*48*24 (EMCal) + 4*48*8 (EMCal/DCal 1/3) + 6*32*24 (DCal)...
Definition: AliAnalysisTaskEMCALClusterize.h:366
AliAnalysisTaskEMCALClusterize::fTCardCorrInduceEnerFracMin
Float_t fTCardCorrInduceEnerFracMin[22]
In case fTCardCorrInduceEnerFracP1 is non null, restrict the minimum fraction of induced energy per S...
Definition: AliAnalysisTaskEMCALClusterize.h:432
AliAnalysisTaskEMCALClusterize::fEMCALNcellsCut
Int_t fEMCALNcellsCut
At least an EMCAL cluster with fNCellsCut cells over fEnergyCut.
Definition: AliAnalysisTaskEMCALClusterize.h:404
AliAnalysisTaskEMCALClusterize
Reclusterize EMCal clusters, put them in a new branch for other following analysis.
Definition: AliAnalysisTaskEMCALClusterize.h:40
AliAnalysisTaskEMCALClusterize::fTCardCorrCellsNew
Bool_t fTCardCorrCellsNew[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells, that before had no signal.
Definition: AliAnalysisTaskEMCALClusterize.h:425
AliAnalysisTaskEMCALClusterize::fLoadGeomMatrices
Bool_t fLoadGeomMatrices
Matrices set from configuration, not get from geometry.root or from ESDs/AODs.
Definition: AliAnalysisTaskEMCALClusterize.h:331
AliAnalysisTaskEMCALClusterize::fEMCALEnergyCut
Float_t fEMCALEnergyCut
At least an EMCAL cluster with this energy in the event.
Definition: AliAnalysisTaskEMCALClusterize.h:403
AliAnalysisTaskEMCALClusterize::fRandomizeTCard
Bool_t fRandomizeTCard
Use random induced energy.
Definition: AliAnalysisTaskEMCALClusterize.h:437
AliAnalysisTaskEMCALClusterize::fUseAliCentrality
Bool_t fUseAliCentrality
Use the centrality estimator from AliCentrality or AliMultSelection.
Definition: AliAnalysisTaskEMCALClusterize.h:399
AliAnalysisTaskEMCALClusterize::fTCardCorrInduceEnerFracWidth
Float_t fTCardCorrInduceEnerFracWidth[4][22]
Induced energy loss gauss witdth on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ ...
Definition: AliAnalysisTaskEMCALClusterize.h:430
AliAnalysisTaskEMCALClusterize::fOrgClusterCellId
Int_t fOrgClusterCellId[fgkNEMCalCells]
Array ID of cluster to wich the cell belongs in unmodified clusters.
Definition: AliAnalysisTaskEMCALClusterize.h:368
AliAnalysisTaskEMCALClusterize::fSelectCellMinFrac
Float_t fSelectCellMinFrac
Min fraction of cell energy after unfolding cut.
Definition: AliAnalysisTaskEMCALClusterize.h:383
AliAnalysisTaskEMCALClusterize::fRemapMCLabelForAODs
Bool_t fRemapMCLabelForAODs
Remap AOD cells MC label. Needed in old AOD productions.
Definition: AliAnalysisTaskEMCALClusterize.h:416