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