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