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