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