AliPhysics  aaf9c62 (aaf9c62)
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, status = 0;
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, status) )
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, status = 0;
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, status))
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  Int_t status = 0;
1206  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1207  {
1208  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1209  fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
1210  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1211 
1213  {
1214  calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
1215  }
1216 
1217  if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi, status))
1218  {
1219  // Channel is not declared as bad
1220  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
1221  eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
1222  }
1223  else
1224  {
1225  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
1226  }
1227  }
1228 
1229  aodEMcells.Sort();
1230 }
1231 
1232 //__________________________________________________
1234 //__________________________________________________
1236 {
1237  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
1238  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
1239 
1240  Double_t pos[3] ;
1241  Double_t covVtx[6];
1242  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
1243 
1244  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
1245  if(!header) AliFatal("Not a standard AOD");
1246  header->SetRunNumber(fRun);
1247 
1248  if(esdevent)
1249  {
1250  TTree* tree = fInputHandler->GetTree();
1251  if (tree)
1252  {
1253  TFile* file = tree->GetCurrentFile();
1254  if (file) header->SetESDFileName(file->GetName());
1255  }
1256  }
1257  else if (aodevent) {
1258  AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
1259  if(!aodheader) AliFatal("Not a standard AOD");
1260  header->SetESDFileName(aodheader->GetESDFileName());
1261  }
1262 
1263  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
1264  header->SetOrbitNumber(fEvent->GetOrbitNumber());
1265  header->SetPeriodNumber(fEvent->GetPeriodNumber());
1266  header->SetEventType(fEvent->GetEventType());
1267 
1268  // Centrality
1269  if(fUseAliCentrality)
1270  {
1271  if(fEvent->GetCentrality())
1272  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
1273  else
1274  header->SetCentrality(0);
1275  }
1276 
1277  // Trigger
1278  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
1279 
1280  header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
1281 
1282  header->SetTriggerMask(fEvent->GetTriggerMask());
1283  header->SetTriggerCluster(fEvent->GetTriggerCluster());
1284 
1285  if (esdevent)
1286  {
1287  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
1288  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
1289  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
1290  }
1291  else if (aodevent)
1292  {
1293  header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
1294  header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
1295  header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
1296  }
1297 
1298  header->SetMagneticField(fEvent->GetMagneticField());
1299  //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
1300 
1301  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
1302  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
1303  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
1304  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
1305  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
1306 
1307  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
1308  Float_t diamcov[3];
1309  fEvent->GetDiamondCovXY(diamcov);
1310  header->SetDiamond(diamxy,diamcov);
1311  if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
1312  else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
1313 
1314  //
1315  Int_t nVertices = 1 ;/* = prim. vtx*/;
1316  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1317 
1318  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1319 
1320  // Access to the AOD container of vertices
1321  TClonesArray &vertices = *(AODEvent()->GetVertices());
1322  Int_t jVertices=0;
1323 
1324  // Add primary vertex. The primary tracks will be defined
1325  // after the loops on the composite objects (V0, cascades, kinks)
1326  fEvent->GetPrimaryVertex()->GetXYZ(pos);
1327  Float_t chi = 0;
1328  if (esdevent)
1329  {
1330  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1331  chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
1332  }
1333  else if (aodevent)
1334  {
1335  aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1336  chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
1337  }
1338 
1339  AliAODVertex * primary = new(vertices[jVertices++])
1340  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
1341  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
1342  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
1343 }
1344 
1345 //___________________________________________________________
1349 //___________________________________________________________
1351 {
1352  // First recalculate track-matching for the new clusters
1353  if(fDoTrackMatching)
1354  {
1356  }
1357 
1358  // Put the new clusters in the AOD list
1359 
1360  Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
1361 
1362  for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
1363  {
1364  AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1365 
1366  newCluster->SetID(i);
1367 
1368  // Correct cluster energy non linearity
1369  newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
1370 
1371  //Add matched track
1372  if(fDoTrackMatching)
1373  {
1374  Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1375  if(trackIndex >= 0)
1376  {
1377  newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1378  AliDebug(2,Form("Matched Track index %d to new cluster %d",trackIndex,i));
1379  }
1380 
1381  Float_t dR = 999., dZ = 999.;
1382  fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dZ,dR);
1383  newCluster->SetTrackDistance(dR,dZ);
1384  }
1385  else
1386  {// Assign previously assigned matched track in reco, very very rough
1387  Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1388  newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1389  }
1390 
1391  //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1392  //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1393  //printf("\n");
1394  //
1395  //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",
1396  // i, newCluster->GetID(), newCluster->E(), newCluster->GetTOF()*1e9, newCluster->GetNCells(), newCluster->GetNLabels(), newCluster->GetLabel());
1397  //
1398  //for(Int_t imc = 0; imc < newCluster->GetNLabels(); imc++)
1399  // printf("%d) Label %d, E dep frac %0.2f; ",
1400  // imc,newCluster->GetLabelAt(imc),newCluster->GetClusterMCEdepFraction(imc));
1401  //if(newCluster->GetNLabels() > 0) printf("\n");
1402 
1403  // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1405 
1406  new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
1407 
1408  AliDebug(2,Form("New cluster %d of %d, energy %f, mc label %d",
1409  newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()));
1410 
1411  } // cluster loop
1412 
1413  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1414 }
1415 
1416 
1417 //________________________________________________________________
1419 //________________________________________________________________
1421 {
1422  if(fUseAliCentrality)
1423  {
1424  if(GetCentrality()) return GetCentrality()->GetCentralityPercentile(fCentralityClass) ;
1425  else return -1. ;
1426  }
1427  else
1428  {
1429  if(GetMultSelCen()) return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass,kTRUE) ;
1430  else return -1. ;
1431  }
1432 
1433  return -1.;
1434 }
1435 
1436 //_______________________________________________
1438 //_______________________________________________
1440 {
1441  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1442  {
1443  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null");
1444  return TString("");
1445  }
1446 
1447  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1448  {
1449  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null");
1450  return TString("");
1451  }
1452 
1453  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1454  if (pass.Contains("ass1")) return TString("pass1");
1455  else if (pass.Contains("ass2")) return TString("pass2");
1456  else if (pass.Contains("ass3")) return TString("pass3");
1457  else if (pass.Contains("ass4")) return TString("pass4");
1458  else if (pass.Contains("ass5")) return TString("pass5");
1459  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1460  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1461  {
1462  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1463  return TString("pass1");
1464  }
1465  else if (pass.Contains("LHC14a1a"))
1466  {
1467  AliInfo("Enable EMCal energy calibration for this MC production!!");
1468 
1470 
1471  return TString("LHC14a1a");
1472  }
1473 
1474  // No condition fullfilled, give a default value
1475  AliInfo("Pass number string not found");
1476 
1477  return TString("");
1478 }
1479 
1480 //_________________________________________
1483 //_________________________________________
1485 {
1486  if(fDebug >=0) (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1487 
1488  fOADBSet = kFALSE;
1489 
1490  fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1491 
1492  if(!fRecParam) fRecParam = new AliEMCALRecParam;
1493  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1494 
1495  if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1496  if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1497  if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1498  fRejectBelowThreshold = kFALSE;
1499 
1500  // Centrality
1501  if(fCentralityClass == "") fCentralityClass = "V0M";
1502 
1503  if (fOCDBpath == "") fOCDBpath = "raw://" ;
1504  if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1505 
1506  if(gROOT->LoadMacro(fConfigName) >=0)
1507  {
1508  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
1509  AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1510  fGeomName = clus->fGeomName;
1512  fOCDBpath = clus->fOCDBpath;
1513  fAccessOCDB = clus->fAccessOCDB;
1514  fRecParam = clus->fRecParam;
1515  fJustUnfold = clus->fJustUnfold;
1516  fFillAODFile = clus->fFillAODFile;
1517  fRecoUtils = clus->fRecoUtils;
1518  fConfigName = clus->fConfigName;
1519  fMaxEvent = clus->fMaxEvent;
1520  fMinEvent = clus->fMinEvent;
1522  fUpdateCell = clus->fUpdateCell;
1524  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1526  fCentralityBin[0] = clus->fCentralityBin[0];
1527  fCentralityBin[1] = clus->fCentralityBin[1];
1529  }
1530 }
1531 
1532 //_______________________________________________________
1535 //_______________________________________________________
1537 {
1538  if (fJustUnfold)
1539  {
1540  // init the unfolding afterburner
1541  delete fUnfolder;
1542  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1543  return;
1544  }
1545 
1546  //First init the clusterizer
1547  delete fClusterizer;
1548  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1549  fClusterizer = new AliEMCALClusterizerv1 (fGeom);
1550  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1551  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1552  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1553  {
1554  fClusterizer = new AliEMCALClusterizerNxN(fGeom);
1555  fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1556  fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1557  }
1558  else
1559  {
1560  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1561  }
1562 
1563  //Now set the parameters
1564  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1565  fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1566  fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1567  fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1568  fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1569  fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1570  fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1571  fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1572  fClusterizer->SetTimeCalibration ( fRecParam->IsTimeCalibrationOn() );
1573  fClusterizer->SetInputCalibrated ( kTRUE );
1574  fClusterizer->SetJustClusters ( kTRUE );
1575 
1576  // Initialize the cluster rec points and digits arrays and get them.
1577  fClusterizer->SetOutput(0);
1578  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1579  fDigitsArr = fClusterizer->GetDigits();
1580 
1581  //In case of unfolding after clusterization is requested, set the corresponding parameters
1582  if(fRecParam->GetUnfold())
1583  {
1584  Int_t i=0;
1585  for (i = 0; i < 8; i++)
1586  {
1587  fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1588  }//end of loop over parameters
1589 
1590  for (i = 0; i < 3; i++)
1591  {
1592  fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1593  fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1594  }//end of loop over parameters
1595 
1596  fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1597  fClusterizer->InitClusterUnfolding();
1598 
1599  }// to unfold
1600 }
1601 
1602 //________________________________________________________________
1606 //________________________________________________________________
1608 {
1609  if(fGeomMatrixSet) return;
1610 
1611  if (!fGeom)
1612  {
1613  if(fGeomName=="")
1614  {
1615  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(fRun);
1616  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fGeom->GetName(),fRun));
1617  }
1618  else
1619  {
1620  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1621  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeomName.Data()));
1622  }
1623 
1624  // Init geometry, I do not like much to do it like this ...
1625  if(fImportGeometryFromFile && !gGeoManager)
1626  {
1627  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1628  {
1629  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1630  if (fRun < 140000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2010.root").data();
1631  else if (fRun < 171000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2011.root").data();
1632  else if (fRun < 198000) fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2012.root").data(); // 2012-2013
1633  else fImportGeometryFilePath = AliDataFile::GetFileNameOADB("EMCAL/geometry_2015.root").data(); // >=2015
1634  }
1635 
1636  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1637 
1638  TGeoManager::Import(fImportGeometryFilePath) ;
1639  }
1640 
1641  AliDebug(1,Form("Init for run=%d",fRun));
1642  if (!gGeoManager) AliDebug(1,"Careful!, gGeoManager not loaded, load misalign matrices");
1643  } // geometry pointer did not exist before
1644 
1645  if(fLoadGeomMatrices)
1646  {
1647  AliInfo("Load user defined EMCAL geometry matrices");
1648 
1649  // OADB if available
1650  AliOADBContainer emcGeoMat("AliEMCALgeo");
1651  if(fOADBFilePath!="")
1652  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1653  else
1654  emcGeoMat.InitFromFile(AliDataFile::GetFileNameOADB("EMCAL/EMCALlocal2master.root").data(),"AliEMCALgeo");
1655 
1656  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(fRun,"EmcalMatrices");
1657 
1658  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1659  {
1660 
1661  if (!fGeomMatrix[mod]) // Get it from OADB
1662  {
1663  AliDebug(2,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
1664  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1665 
1666  fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1667  }
1668 
1669  if(fGeomMatrix[mod])
1670  {
1671  if(DebugLevel() > 1)
1672  fGeomMatrix[mod]->Print();
1673 
1674  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1675  }
1676  else if(gGeoManager)
1677  {
1678  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
1679  fGeom->SetMisalMatrix(fGeom->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
1680  }
1681  else
1682  {
1683  AliError(Form("Alignment atrix for SM %d is not available",mod));
1684  }
1685 
1686  fGeomMatrixSet=kTRUE;
1687 
1688  }//SM loop
1689  }//Load matrices
1690  else if(!gGeoManager)
1691  {
1692  AliInfo("Get geo matrices from data");
1693  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1694  if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1695  {
1696  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
1697  }//AOD
1698  else
1699  {
1700  AliDebug(1,"Load Misaligned matrices");
1701 
1702  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1703 
1704  if(!esd)
1705  {
1706  AliError("This event does not contain ESDs?");
1707  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1708  return;
1709  }
1710 
1711  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1712  {
1713  if(DebugLevel() > 1)
1714  esd->GetEMCALMatrix(mod)->Print();
1715 
1716  if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1717 
1718  }
1719 
1720  fGeomMatrixSet=kTRUE;
1721 
1722  } // ESD
1723  } // Load matrices from Data
1724 }
1725 
1726 //____________________________________________________
1731 //____________________________________________________
1733 {
1734  if(!fRemoveExoticEvents) return kFALSE;
1735 
1736  // Loop on cells
1737 
1738  Float_t totCellE = 0;
1739  Int_t bc = InputEvent()->GetBunchCrossNumber();
1740 
1741  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1742  {
1743  Float_t ecell = 0 ;
1744  Double_t tcell = 0 ;
1745 
1746  Int_t absID = fCaloCells->GetCellNumber(icell);
1747  Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,fCaloCells);
1748  tcell-=fConstantTimeShift*1e-9;// Only for MC simulations done before 2015
1749 
1750  if(accept && !fRecoUtils->IsExoticCell(absID,fCaloCells,bc)) totCellE += ecell;
1751  }
1752 
1753  // TString triggerclasses = event->GetFiredTriggerClasses();
1754  // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1755  // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1756  // return;
1757  //
1758 
1759  //printf("TotE cell %f\n",totCellE);
1760  if(totCellE < 1) return kTRUE;
1761 
1762  return kFALSE;
1763 }
1764 
1765 //________________________________________________________________
1767 //________________________________________________________________
1769 {
1770  if(!fRemoveLEDEvents) return kFALSE;
1771 
1772  //check events of LHC11a period
1773  if(run < 146858 || run > 146860) return kFALSE ;
1774 
1775  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1776  Int_t ncellsSM3 = 0;
1777  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1778  {
1779  if ( fCaloCells->GetAmplitude (icell) > 0.1 &&
1780  fCaloCells->GetCellNumber(icell)/(24*48)==3 ) ncellsSM3++;
1781  }
1782 
1783  TString triggerclasses = fEvent->GetFiredTriggerClasses();
1784 
1785  Int_t ncellcut = 21;
1786  if(triggerclasses.Contains("EMC")) ncellcut = 35;
1787 
1788  if( ncellsSM3 >= ncellcut)
1789  {
1790  AliInfo(Form("Reject event %d with ncells in SM3 %d",(Int_t)Entry(),ncellsSM3));
1791  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1792  return kTRUE;
1793  }
1794 
1795  return kFALSE;
1796 }
1797 
1798 //_______________________________________________________
1807 //_______________________________________________________
1809 (Int_t absId, Int_t absIdRef, Int_t sm, Float_t ampRef, Int_t cellCase)
1810 {
1811  // Check that the cell exists
1812  if( !AcceptCell(absId,0) ) return ;
1813 
1814  // Get the fraction
1815  Float_t frac = fTCardCorrInduceEnerFrac[cellCase][sm] + ampRef * fTCardCorrInduceEnerFracP1[cellCase][sm];
1816 
1817  // Use an absolute minimum and maximum fraction if calculated one is out of range
1818  if ( frac < fTCardCorrInduceEnerFracMin[sm] ) frac = fTCardCorrInduceEnerFracMin[sm];
1819  if ( frac > fTCardCorrInduceEnerFracMax[sm] ) frac = fTCardCorrInduceEnerFracMax[sm];
1820 
1821  AliDebug(1,Form("\t fraction %2.3f",frac));
1822 
1823  // Randomize the induced fraction, if requested
1824  if ( fRandomizeTCard )
1825  {
1826  frac = fRandom.Gaus(frac, fTCardCorrInduceEnerFracWidth[cellCase][sm]);
1827 
1828  AliDebug(1,Form("\t randomized fraction %2.3f",frac));
1829  }
1830 
1831  // If too small or negative, do nothing else
1832  if ( frac < 0.0001 ) return;
1833 
1834  // Calculate induced energy
1835  Float_t inducedE = fTCardCorrInduceEner[cellCase][sm] + ampRef * frac;
1836 
1837  // Check if we induce too much energy, in such case use a constant value
1838  if ( fTCardCorrMaxInduced < inducedE ) inducedE = fTCardCorrMaxInduced;
1839 
1840  AliDebug(1,Form("\t induced E %2.3f",inducedE));
1841 
1842  // Add the induced energy, check if cell existed
1843  // Check that the induced+amp is large enough to avoid extra linearity effects
1844  // typically of the order of the clusterization cell energy cut
1845  // But if it is below 1 ADC, typically 10 MeV, also do it, to match Beam test linearity
1846  Float_t amp = fCaloCells->GetCellAmplitude(absId) ;
1847  if ( (amp+inducedE) > fTCardCorrMinInduced || inducedE < fTCardCorrMaxInducedLowE )
1848  {
1849  fTCardCorrCellsEner[absId] += inducedE;
1850 
1851  // If original energy of cell was null, create new one
1852  if ( amp < 0.01 ) fTCardCorrCellsNew[absId] = kTRUE;
1853  }
1854  else return ;
1855 
1856  // Subtract the added energy to main cell, if energy conservation is requested
1858  fTCardCorrCellsEner[absIdRef] -= inducedE;
1859 }
1860 
1861 //_______________________________________________________
1864 //_______________________________________________________
1866 {
1867  Int_t id = -1;
1868  Float_t amp = -1;
1869 
1870  // Loop on all cells with signal
1871  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1872  {
1873  id = fCaloCells->GetCellNumber(icell);
1874  amp = fCaloCells->GetAmplitude (icell); // fCaloCells->GetCellAmplitude(id);
1875 
1876  if ( amp <= fTCardCorrMinAmp ) continue ;
1877 
1878  //
1879  // First get the SM, col-row of this tower
1880  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1881  fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
1882  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1883 
1884  //
1885  // Determine randomly if we want to create a correlation for this cell,
1886  // depending the SM number of the cell
1887  if ( fTCardCorrInduceEnerProb[imod] < 1 )
1888  {
1889  Float_t rand = fRandom.Uniform(0, 1);
1890 
1891  if ( rand > fTCardCorrInduceEnerProb[imod] )
1892  {
1893  AliDebug(1,Form("Do not difuse E of cell %d, sm %d, amp %2.2f: SM fraction %2.2f > %2.2f",
1894  id,imod,amp,fTCardCorrInduceEnerProb[imod],rand));
1895  continue;
1896  }
1897  }
1898 
1899  AliDebug(1,Form("Reference cell absId %d, iEta %d, iPhi %d, sm %d, amp %2.2f",id,ieta,iphi,imod,amp));
1900 
1901  //
1902  // Get the absId of the cells in the cross and same T-Card
1903  Int_t absIDup = -1;
1904  Int_t absIDdo = -1;
1905  Int_t absIDlr = -1;
1906  Int_t absIDuplr = -1;
1907  Int_t absIDdolr = -1;
1908  Int_t absIDup2 = -1;
1909  Int_t absIDup2lr = -1;
1910  Int_t absIDdo2 = -1;
1911  Int_t absIDdo2lr = -1;
1912 
1913  // Only 2 columns in the T-Card, +1 for even and -1 for odd with respect reference cell
1914  Int_t colShift = 0;
1915  if ( (ieta%2) && ieta <= AliEMCALGeoParams::fgkEMCALCols-1 ) colShift = -1;
1916  if ( !(ieta%2) && ieta >= 0 ) colShift = +1;
1917 
1918  absIDlr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+colShift);
1919 
1920  // Check up / down cells from reference cell not out of SM and in same T-Card
1921  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1 )
1922  {
1923  absIDup = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
1924  absIDuplr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta+colShift);
1925  }
1926 
1927  if ( iphi > 0 )
1928  {
1929  absIDdo = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
1930  absIDdolr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta+colShift);
1931  }
1932 
1933  // Check 2 up / 2 down cells from reference cell not out of SM
1934  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-2 )
1935  {
1936  absIDup2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta);
1937  absIDup2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta+colShift);
1938  }
1939 
1940  if ( iphi > 1 )
1941  {
1942  absIDdo2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta);
1943  absIDdo2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta+colShift);
1944  }
1945 
1946  // In same T-Card?
1947  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+1)/8) ) { absIDup = -1 ; absIDuplr = -1 ; }
1948  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-1)/8) ) { absIDdo = -1 ; absIDdolr = -1 ; }
1949  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+2)/8) ) { absIDup2 = -1 ; absIDup2lr = -1 ; }
1950  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-2)/8) ) { absIDdo2 = -1 ; absIDdo2lr = -1 ; }
1951 
1952  // Calculate induced energy to T-Card cells
1953 
1954  AliDebug(1,Form("cell up %d:" ,absIDup));
1955  CalculateInducedEnergyInTCardCell(absIDup , id, imod, amp, 0);
1956  AliDebug(1,Form("cell down %d:",absIDdo));
1957  CalculateInducedEnergyInTCardCell(absIDdo , id, imod, amp, 0);
1958 
1959  AliDebug(1,Form("cell up left-right %d:" ,absIDuplr));
1960  CalculateInducedEnergyInTCardCell(absIDuplr , id, imod, amp, 1);
1961  AliDebug(1,Form("cell down left-right %d:",absIDdolr));
1962  CalculateInducedEnergyInTCardCell(absIDdolr , id, imod, amp, 1);
1963 
1964  AliDebug(1,Form("cell left-right %d:",absIDlr));
1965  CalculateInducedEnergyInTCardCell(absIDlr , id, imod, amp, 2);
1966 
1967  AliDebug(1,Form("cell up 2nd row %d:" ,absIDup2));
1968  CalculateInducedEnergyInTCardCell(absIDup2 , id, imod, amp, 3);
1969  AliDebug(1,Form("cell down 2nd row %d:",absIDdo2));
1970  CalculateInducedEnergyInTCardCell(absIDdo2 , id, imod, amp, 3);
1971 
1972  AliDebug(1,Form("cell up left-right 2nd row %d:" ,absIDup2lr));
1973  CalculateInducedEnergyInTCardCell(absIDup2lr, id, imod, amp, 3);
1974  AliDebug(1,Form("cell down left-right 2nd row %d:",absIDdo2lr));
1975  CalculateInducedEnergyInTCardCell(absIDdo2lr, id, imod, amp, 3);
1976 
1977  } // cell loop
1978 
1979 }
1980 
1981 //_______________________________________________________
1983 //_______________________________________________________
1985 {
1986  AliInfo(Form("Geometry: name <%s>, matrix set <%d>, load matrix <%d>, import geo <%d> from path <%s>",
1988 
1989  if ( fAccessOCDB ) AliInfo(Form("OCDB path name <%s>", fOCDBpath.Data()));
1990  if ( fAccessOADB ) AliInfo(Form("OADB path name <%s>", fOADBFilePath.Data()));
1991 
1992  if ( fInputCaloCellsName.Length() > 0 )
1993  AliInfo(Form("Input CaloCells <%s>", fInputCaloCellsName.Data()));
1994 
1995  AliInfo(Form("Just Unfold clusters <%d>, new clusters list name <%s>, new cells name <%s>",
1997 
1998  if ( fFillAODFile ) AliInfo(Form("Fill new AOD file with: header <%d>, cells <%d>",fFillAODHeader,fFillAODCaloCells));
1999 
2000  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>",
2002 
2003  AliInfo(Form("Reject events out of range: %d < N event < %d, LED <%d>, exotics <%d>",
2005 
2006  if (fCentralityBin[0] != -1 && fCentralityBin[1] != -1 )
2007  AliInfo(Form("Centrality bin [%2.2f,%2.2f], class <%s>, use AliCentrality? <%d>",
2009 
2010  if ( fSelectEMCALEvent )
2011  AliInfo(Form("Select events with signal in EMCal: E min <%2.2f>, n cell min <%d>", fEMCALEnergyCut, fEMCALNcellsCut));
2012 
2013  AliInfo(Form("MC label from cluster <%d>, Use EdepFrac <%d>, remap AODs <%d>",
2015 }
2016 
2017 //_______________________________________________________
2019 //_______________________________________________________
2021 {
2022  if(!fTCardCorrEmulation)
2023  {
2024  AliInfo("T-Card emulation not activated");
2025  return;
2026  }
2027 
2028  AliInfo(Form("T-Card emulation activated, energy conservation <%d>, randomize E <%d>, induced energy parameters:",
2030 
2031  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",
2033 
2034  for(Int_t ism = 0; ism < 22; ism++)
2035  {
2036  printf("\t sm %d, fraction %2.3f, E frac abs min %2.3e max %2.3e \n",
2038 
2039  for(Int_t icell = 0; icell < 4; icell++)
2040  {
2041  printf("\t \t cell type %d, c %2.4e, p0 %2.4e, p1 %2.4e, sigma %2.4e \n",
2042  icell,fTCardCorrInduceEner[icell][ism],fTCardCorrInduceEnerFrac[icell][ism],
2044  }
2045  }
2046 }
2047 
2048 //_______________________________________________________
2052 //_______________________________________________________
2054 {
2055  Int_t j = 0;
2056  for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
2057  {
2058  AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
2059 
2060  const Int_t ncells = recPoint->GetMultiplicity();
2061  Int_t ncellsTrue = 0;
2062 
2063  if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
2064 
2065  // cells and their amplitude fractions
2066  UShort_t absIds[ncells];
2067  Double32_t ratios[ncells];
2068 
2069  //For later check embedding
2070  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2071 
2072  Float_t clusterE = 0;
2073  for (Int_t c = 0; c < ncells; c++)
2074  {
2075  AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
2076 
2077  absIds[ncellsTrue] = digit->GetId();
2078  ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
2079 
2080  if ( !fRecParam->GetUnfold() && (ratios[ncellsTrue] > 1 || ratios[ncellsTrue] < 1) )
2081  AliWarning(Form("recpoint cell E %2.3f but digit E %2.3f and no unfolding", recPoint->GetEnergiesList()[c], digit->GetAmplitude()));
2082 
2083  // In case of unfolding, remove digits with unfolded energy too low
2084  if(fSelectCell)
2085  {
2086  if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
2087  {
2088 
2089  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
2090  recPoint->GetEnergiesList()[c],digit->GetAmplitude()));
2091  continue;
2092  } // if cuts
2093  }// Select cells
2094 
2095  //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
2096  clusterE +=recPoint->GetEnergiesList()[c];
2097 
2098  // In case of embedding, fill ratio with amount of signal,
2099  if (aodIH && aodIH->GetMergeEvents())
2100  {
2101  //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
2102  AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
2103  AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
2104 
2105  Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2106  //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2107  Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
2108  //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
2109 
2110  if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
2111  //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
2112 
2113  }//Embedding
2114 
2115  ncellsTrue++;
2116 
2117  }// cluster cell loop
2118 
2119  if (ncellsTrue < 1)
2120  {
2121  AliDebug(2,Form("Skipping cluster with no cells avobe threshold E = %f, ncells %d",
2122  recPoint->GetEnergy(), ncells));
2123  continue;
2124  }
2125 
2126  //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
2127 
2128  if(clusterE < fRecParam->GetClusteringThreshold())
2129  {
2130  AliDebug(2,Form("Remove cluster with energy below seed threshold %f",clusterE));
2131  continue;
2132  }
2133 
2134  TVector3 gpos;
2135  Float_t g[3];
2136 
2137  // calculate new cluster position
2138 
2139  recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
2140  recPoint->GetGlobalPosition(gpos);
2141  gpos.GetXYZ(g);
2142 
2143  // create a new cluster
2144 
2145  (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
2146  AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
2147  j++;
2148  clus->SetType(AliVCluster::kEMCALClusterv1);
2149  clus->SetE(clusterE);
2150  clus->SetPosition(g);
2151  clus->SetNCells(ncellsTrue);
2152  clus->SetCellsAbsId(absIds);
2153  clus->SetCellsAmplitudeFraction(ratios);
2154  clus->SetChi2(-1); //not yet implemented
2155  clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
2156  clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
2157  clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
2158 
2159  if(ncells == ncellsTrue)
2160  {
2161  Float_t elipAxis[2];
2162  recPoint->GetElipsAxis(elipAxis);
2163  clus->SetM02(elipAxis[0]*elipAxis[0]) ;
2164  clus->SetM20(elipAxis[1]*elipAxis[1]) ;
2165  clus->SetDispersion(recPoint->GetDispersion());
2166  }
2167  else if(fSelectCell)
2168  {
2169  // In case some cells rejected, in unfolding case, recalculate
2170  // shower shape parameters and position
2171  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
2172 
2173  AliVCaloCells* cells = 0x0;
2174  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
2175  else cells = InputEvent()->GetEMCALCells();
2176 
2180  }
2181 
2182  // MC
2183 
2186  else
2187  {
2188  //
2189  // Normal case, trust what the clusterizer has found
2190  //
2191  Int_t parentMult = 0;
2192  Int_t *parentList = recPoint->GetParents(parentMult);
2193  Float_t *parentListDE = recPoint->GetParentsDE(); // deposited energy
2194 
2195  clus->SetLabel(parentList, parentMult);
2197  clus->SetClusterMCEdepFractionFromEdepArray(parentListDE);
2198 
2199  //
2200  // Set the cell energy deposition fraction map:
2201  //
2202  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
2203  {
2204  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
2205 
2206  // Get the digit that originated this cell cluster
2207 // AliVCaloCells* cells = 0x0;
2208 // if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
2209 // else cells = InputEvent()->GetEMCALCells();
2210 
2211  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
2212  {
2213  Int_t idigit = fCellLabels[absIds[icell]];
2214 
2215  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
2216 
2217  // Find the 4 MC labels that contributed to the cluster and their
2218  // deposited energy in the current digit
2219 
2220  mcEdepFracPerCell[icell] = 0; // init
2221 
2222  Int_t nparents = dig->GetNiparent();
2223  if ( nparents > 0 )
2224  {
2225  Int_t digLabel =-1 ;
2226  Float_t edep = 0 ;
2227  Float_t edepTot = 0 ;
2228  Float_t mcEDepFrac[4] = {0,0,0,0};
2229 
2230  // all parents in digit
2231  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2232  {
2233  digLabel = dig->GetIparent (jndex+1);
2234  edep = dig->GetDEParent(jndex+1);
2235  edepTot += edep;
2236 
2237  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
2238  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
2239  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
2240  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
2241  } // all prarents in digit
2242 
2243  // Divide energy deposit by total deposited energy
2244  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
2245  if(edepTot > 0.01)
2246  {
2247  mcEdepFracPerCell[icell] = clus->PackMCEdepFraction(mcEDepFrac);
2248  }
2249  } // at least one parent label in digit
2250  } // cell in cluster loop
2251 
2252  clus->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
2253 
2254  delete [] mcEdepFracPerCell;
2255 
2256  } // at least one parent in cluster, do the cell primary packing
2257  }
2258  } // recPoints loop
2259 }
2260 
2261 //_____________________________________________________________________
2264 //_____________________________________________________________________
2266 {
2267  if(label < 0) return ;
2268 
2269  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
2270  if(!evt) return ;
2271 
2272  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
2273  if(!arr) return ;
2274 
2275  if(label < arr->GetEntriesFast())
2276  {
2277  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
2278  if(!particle) return ;
2279 
2280  if(label == particle->Label()) return ; // label already OK
2281  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
2282  }
2283  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
2284 
2285  // loop on the particles list and check if there is one with the same label
2286  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
2287  {
2288  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
2289  if(!particle) continue ;
2290 
2291  if(label == particle->Label())
2292  {
2293  label = ind;
2294  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
2295  return;
2296  }
2297  }
2298 
2299  label = -1;
2300 
2301  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
2302 }
2303 
2304 //________________________________________________
2306 //________________________________________________
2308 {
2309  for(Int_t j = 0; j < fgkNEMCalCells; j++)
2310  {
2311  fOrgClusterCellId[j] = -1;
2312  fCellLabels[j] = -1;
2313  fCellSecondLabels[j] = -1;
2314  fCellTime[j] = 0.;
2315  fCellMatchdEta[j] = -999;
2316  fCellMatchdPhi[j] = -999;
2317 
2318  fTCardCorrCellsEner[j] = 0.;
2319  fTCardCorrCellsNew [j] = kFALSE;
2320  }
2321 }
2322 
2323 //_____________________________________________________________________________________________________
2327 //_____________________________________________________________________________________________________
2329  AliAODCaloCluster * clus)
2330 {
2331  Int_t parentMult = 0;
2332  Int_t *parentList = recPoint->GetParents(parentMult);
2333  clus->SetLabel(parentList, parentMult);
2334 
2335  //Write the second major contributor to each MC cluster.
2336  Int_t iNewLabel ;
2337  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2338  {
2339 
2340  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2341  if(idCell>=0)
2342  {
2343  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
2344  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
2345  {
2346  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
2347  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
2348  }
2349  if (iNewLabel == 1)
2350  {
2351  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
2352  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
2353  {
2354  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
2355  }
2356 
2357  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
2358  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
2359  delete [] newLabelArray;
2360  }
2361  } // positive cell number
2362  }
2363 }
2364 
2365 //___________________________________________________________________________________________________
2369 //___________________________________________________________________________________________________
2371 {
2372  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
2373  clArray.Reset();
2374  Int_t nClu = 0;
2375  Int_t nLabTotOrg = 0;
2376  Float_t emax = -1;
2377  Int_t idMax = -1;
2378 
2379  AliVEvent * event = fEvent;
2380 
2381  // In case of embedding the MC signal is in the added event
2382  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2383  if(aodIH && aodIH->GetEventToMerge()) //Embedding
2384  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
2385 
2386 
2387  // Find the clusters that originally had the cells
2388  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2389  {
2390  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2391 
2392  if(idCell>=0)
2393  {
2394  Int_t idCluster = fOrgClusterCellId[idCell];
2395 
2396  Bool_t set = kTRUE;
2397  for(Int_t icl =0; icl < nClu; icl++)
2398  {
2399  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
2400  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
2401  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
2402  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
2403  }
2404  if( set && idCluster >= 0)
2405  {
2406  clArray.SetAt(idCluster,nClu++);
2407  //printf("******** idCluster %d \n",idCluster);
2408  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
2409 
2410  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
2411 
2412  //Search highest E cluster
2413  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2414  //printf("\t E %f\n",clOrg->E());
2415  if(emax < clOrg->E())
2416  {
2417  emax = clOrg->E();
2418  idMax = idCluster;
2419  }
2420  }
2421  }
2422  }// cell loop
2423 
2424  if(nClu==0 || nLabTotOrg == 0)
2425  {
2426  //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());
2427  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
2428  //printf("\n");
2429  }
2430 
2431  // Put the first in the list the cluster with highest energy
2432  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
2433  {
2434  Int_t maxIndex = -1;
2435  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
2436  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2437  {
2438  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
2439  }
2440 
2441  if(firstCluster >=0 && idMax >=0)
2442  {
2443  clArray.SetAt(idMax,0);
2444  clArray.SetAt(firstCluster,maxIndex);
2445  }
2446  }
2447 
2448  // Get the labels list in the original clusters, assign all to the new cluster
2449  TArrayI clMCArray(nLabTotOrg) ;
2450  clMCArray.Reset();
2451 
2452  Int_t nLabTot = 0;
2453  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2454  {
2455  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
2456  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
2457  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2458  Int_t nLab = clOrg->GetNLabels();
2459 
2460  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
2461  {
2462  Int_t lab = clOrg->GetLabelAt(iLab) ;
2463  if(lab>=0)
2464  {
2465  Bool_t set = kTRUE;
2466  //printf("\t \t Set Label %d \n", lab);
2467  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
2468  {
2469  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
2470  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
2471  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
2472  }
2473  if( set ) clMCArray.SetAt(lab,nLabTot++);
2474  }
2475  }
2476  }// cluster loop
2477 
2478  // Set the final list of labels
2479 
2480  clus->SetLabel(clMCArray.GetArray(), nLabTot);
2481 
2482 // printf("Final list of labels for new cluster : \n");
2483 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
2484 // {
2485 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
2486 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
2487 // printf(" org %d ",label);
2488 // RemapMCLabelForAODs(label);
2489 // printf(" new %d \n",label);
2490 // }
2491 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
2492 }
2493 
2494 //____________________________________________________________
2495 // Init geometry, create list of output clusters and cells.
2496 //____________________________________________________________
2498 {
2499  // Clusters
2500  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2501 
2502  if(fOutputAODBranchName.Length()==0)
2503  {
2504  fOutputAODBranchName = "newEMCALClustersArray";
2505  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
2506  }
2507 
2509 
2510  if( fFillAODFile )
2511  {
2512  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2513 
2514  //fOutputAODBranch->SetOwner(kFALSE);
2515 
2516  AddAODBranch("TClonesArray", &fOutputAODBranch);
2517  }
2518 
2519  // Cells
2520  if ( fOutputAODBranchName.Length() > 0 )
2521  {
2522  fOutputAODCells = new AliAODCaloCells(fOutputAODCellsName,fOutputAODCellsName,AliAODCaloCells::kEMCALCell);
2523 
2524  if( fFillAODFile ) AddAODBranch("AliAODCaloCells", &fOutputAODCells);
2525  }
2526 }
2527 
2528 //_______________________________________________________________________
2531 //________________________________________________________________________
2533 {
2534  // Update cells only in case re-calibration was done
2535  // or bad map applied or additional T-Card cells added.
2541  !fTCardCorrEmulation ) return;
2542 
2543  const Int_t ncells = fCaloCells->GetNumberOfCells();
2544  const Int_t ndigis = fDigitsArr->GetEntries();
2545 
2546  if ( fOutputAODCellsName.Length() > 0 )
2547  {
2548  fOutputAODCells->DeleteContainer();
2549  fOutputAODCells->CreateContainer(ndigis);
2550  }
2551  else if ( ncells != ndigis ) // update case
2552  {
2553  fCaloCells->DeleteContainer();
2554  fCaloCells->CreateContainer(ndigis);
2555  }
2556 
2557  for (Int_t idigit = 0; idigit < ndigis; ++idigit)
2558  {
2559  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
2560 
2561  Double_t cellAmplitude = digit->GetAmplitude();
2562  Short_t cellNumber = digit->GetId();
2563  Double_t cellTime = digit->GetTime();
2564 
2565  Bool_t highGain = kFALSE;
2566  if( digit->GetType() == AliEMCALDigit::kHG ) highGain = kTRUE;
2567 
2568  // Only for MC
2569  // Get the label of the primary particle that generated the cell
2570  // Assign the particle that deposited more energy
2571  Int_t nparents = digit->GetNiparent();
2572  Int_t cellMcEDepFrac =-1 ;
2573  Float_t cellMcLabel =-1.;
2574  if ( nparents > 0 )
2575  {
2576  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2577  {
2578  if(cellMcEDepFrac >= digit->GetDEParent(jndex+1)) continue ;
2579 
2580  cellMcLabel = digit->GetIparent (jndex+1);
2581  cellMcEDepFrac= digit->GetDEParent(jndex+1);
2582  } // all primaries in digit
2583  } // select primary label
2584 
2585  if ( cellMcEDepFrac < 0 ) cellMcEDepFrac = 0.;
2586 
2587  if ( fUpdateCell )
2588  fCaloCells ->SetCell(idigit, cellNumber, cellAmplitude, cellTime, cellMcLabel, cellMcEDepFrac, highGain);
2589  else
2590  fOutputAODCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime, cellMcLabel, cellMcEDepFrac, highGain);
2591  }
2592 
2593  if ( ncells != ndigis )
2594  {
2595  if ( fUpdateCell )
2596  fCaloCells ->Sort();
2597  else
2598  fOutputAODCells->Sort();
2599  }
2600 }
2601 
2602 //_______________________________________________________
2613 //_______________________________________________________
2615 {
2616  //-------
2617  // Step 1
2618 
2619  // Remove the contents of AOD branch output list set in the previous event
2620  fOutputAODBranch->Clear("C");
2621 
2622  LoadBranches();
2623 
2624  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
2625  // If we need a centrality bin, we select only those events in the corresponding bin.
2626  if( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
2627  {
2628  Float_t cen = GetEventCentrality();
2629  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
2630  }
2631 
2632  // Intermediate array with new clusters : init the array only once or clear from previous event
2633  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
2634  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
2635 
2636 
2637  // In case of analysis in multiple runs, check the OADB again
2638  if ( InputEvent()->GetRunNumber() != fRun )
2639  {
2640  fRun = InputEvent()->GetRunNumber();
2641 
2642  fOADBSet = kFALSE; // recover the OADB for this run
2643 
2644  AliInfo(Form("Set run to %d",fRun));
2645  }
2646 
2647  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
2648 
2649  // Get the event, do some checks and settings
2650  CheckAndGetEvent() ;
2651 
2652  if (!fEvent)
2653  {
2654  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
2655  return ;
2656  }
2657 
2658  // Init pointers, geometry, clusterizer, ocdb, aodb
2659 
2660  if(fAccessOCDB) AccessOCDB();
2661  if(fAccessOADB) AccessOADB(); // only once
2662 
2664 
2665  // Print once the analysis parameters
2666  if ( fDebug > 0 || !fPrintOnce )
2667  {
2668  //fRecParam->Print("reco"); // AliInfo not printed ...
2669 
2670  PrintParam();
2671 
2672  PrintTCardParam();
2673 
2674  fPrintOnce = kTRUE;
2675  }
2676 
2677  //-------
2678  // Step 2
2679 
2680  // Make clusters
2682  else ClusterizeCells() ;
2683 
2684  //-------
2685  // Step 3
2686 
2688 
2689  if ( fUpdateCell || fOutputAODCellsName.Length() > 0 )
2690  UpdateCells();
2691 }
2692 
2693 
Bool_t IsRunDepRecalibrationOn() const
void RecalculateCellLabelsRemoveAddedGenerator(Int_t absID, AliVCluster *clus, AliMCEvent *mc, Float_t &amp, TArrayI &labeArr, TArrayF &eDepArr) const
TString fOutputAODCellsName
New of output cells AOD branch name.
TH1F * GetEMCALChannelTimeRecalibrationFactors(Int_t bc) const
void FillAODCaloCells()
Put calo cells in standard branch.
Float_t fTCardCorrMinAmp
Minimum cell energy to induce signal on adjacent cells.
Bool_t IsL1PhaseInTimeRecalibrationOn() const
AliVCaloCells * fCaloCells
! CaloCells container
Bool_t fTCardCorrEmulation
Activate T-Card cells energy correlation.
Float_t fTCardCorrInduceEnerFracP1[4][22]
Induced energy loss gauss fraction param1 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
double Double_t
Definition: External.C:58
virtual void ResetArrays()
Reset arrays containing information for all possible cells.
TH2F * GetEMCALChannelRecalibrationFactors(Int_t iSM) const
Definition: External.C:236
TString GetPass()
Get or guess pass number/string from path of filename.
Bool_t fFillAODHeader
Copy header to standard branch.
TClonesArray * fDigitsArr
! Digits array
Bool_t fTCardCorrClusEnerConserv
When making correlation, subtract from the reference cell the induced energy on the neighbour cells...
void InitEMCALL1PhaseInTimeRecalibration()
AliEMCALClusterizer * fClusterizer
! EMCAL clusterizer
TH1C * GetEMCALL1PhaseInTimeRecalibrationForAllSM() const
void SetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow, Double_t c=1)
Bool_t fJustUnfold
Just unfold, do not recluster.
Bool_t fImportGeometryFromFile
Import geometry settings in geometry.root file.
Bool_t fSelectCell
Reject cells from cluster if energy is too low and recalculate position/energy and other...
TString fConfigName
Name of analysis configuration file.
Float_t fTCardCorrInduceEnerFracMax[22]
In case fTCardCorrInduceEnerFracP1 is non null, restrict the maximum fraction of induced energy per S...
TCanvas * c
Definition: TestFitELoss.C:172
void PrintTCardParam()
Print parameters for T-Card correlation emulation.
Int_t fCellSecondLabels[fgkNEMCalCells]
Array with Second MC label to be passed to digit.
void SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint *recPoint, AliAODCaloCluster *clus)
Bool_t IsBadChannelsRemovalSwitchedOn() const
Bool_t IsExoticCell(Int_t absId, AliVCaloCells *cells, Int_t bc=-1)
void SetNonLinearityFunction(Int_t fun)
Float_t fCellMatchdEta[fgkNEMCalCells]
Array with cluster-track dPhi.
TClonesArray * fOutputAODBranch
! AOD Branch with output clusters
TString fOADBFilePath
Default path $ALICE_PHYSICS/OADB/EMCAL, if needed change.
Float_t fTCardCorrCellsEner[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells.
Bool_t IsRecalibrationOn() const
Some utilities for cluster and cell treatment.
TString fImportGeometryFilePath
path fo geometry.root file
Float_t fTCardCorrInduceEnerFrac[4][22]
Induced energy loss gauss fraction param0 on 0-same row, diff col, 1-up/down cells left/right col 2-l...
TObjArray * GetEMCALL1PhaseInTimeRecalibrationArray() const
void ConfigureEMCALRecoUtils(Bool_t bMC=kFALSE, Bool_t bExotic=kTRUE, Bool_t bNonLin=kFALSE, Bool_t bRecalE=kTRUE, Bool_t bBad=kTRUE, Bool_t bRecalT=kTRUE, Int_t debug=-1)
Bool_t fPrintOnce
Print once analysis parameters.
Bool_t fSelectEMCALEvent
Process the event if there is some high energy cluster.
Float_t fTCardCorrMaxInduced
Maximum induced energy signal on adjacent cells.
Float_t GetEventCentrality() const
Get centrality/multiplicity percentile.
void RecalculateClusterDistanceToBadChannel(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
TH2I * GetEMCALChannelStatusMap(Int_t iSM) const
Float_t GetEMCALChannelRecalibrationFactor(Int_t iSM, Int_t iCol, Int_t iRow) const
Float_t fTCardCorrMaxInducedLowE
Maximum value of induced energy signal that is always accepted, order of ADC, tipically 10 MeV...
Bool_t fInputFromFilter
Get the input from AODs from the filter.
TString fOutputAODBranchName
New of output clusters AOD branch.
Int_t GetMatchedTrackIndex(Int_t clsIndex)
AliEMCALGeometry * fGeom
EMCAL geometry.
void SwitchOnL1PhaseInTimeRecalibration()
void ClusterUnfolding()
Take the event clusters and unfold them.
int Int_t
Definition: External.C:63
Bool_t fAccessOADB
Get calibration from OADB for EMCAL.
Bool_t fRemoveLEDEvents
Remove LED events, use only for LHC11a.
Bool_t AcceptCalibrateCell(Int_t absId, Int_t bc, Float_t &amp, Double_t &time, AliVCaloCells *cells)
unsigned int UInt_t
Definition: External.C:33
Bool_t GetEMCALChannelStatus(Int_t iSM, Int_t iCol, Int_t iRow, Int_t &status) const
AliEMCALRecoUtils * fRecoUtils
Access to factorized reconstruction algorithms.
void RecalculateClusterShowerShapeParameters(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *cluster)
float Float_t
Definition: External.C:68
Bool_t fRejectBelowThreshold
split (false-default) or reject (true) cell energy below threshold after UF
Float_t fConstantTimeShift
Apply a 600 ns time shift in case of simulation, shift in ns.
void SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster *clus)
Bool_t fUpdateCell
On/Off the upate of the CaloCells container.
Float_t fSelectCellMinE
Min energy cell threshold, after unfolding.
Bool_t fOADBSet
AODB parameters already set.
AliAODCaloCells * fOutputAODCells
! AOD Branch with output cells
Bool_t IsTimeRecalibrationOn() const
Float_t fTCardCorrInduceEner[4][22]
Induced energy loss gauss constant on 0-same row, diff col, 1-up/down cells left/right col 2-left/rig...
Int_t fCellLabels[fgkNEMCalCells]
Array with MC label to be passed to digit.
Bool_t fRemoveExoticEvents
Remove exotic events.
Bool_t fDoTrackMatching
On/Off the matching recalculation to speed up analysis in PbPb.
Float_t fCellMatchdPhi[fgkNEMCalCells]
Array with cluster-track dEta.
Float_t fTCardCorrInduceEnerProb[22]
Probability to induce energy loss per SM.
void RecalibrateClusterEnergy(const AliEMCALGeometry *geom, AliVCluster *cluster, AliVCaloCells *cells, Int_t bc=-1)
void GetMatchedResiduals(Int_t clsIndex, Float_t &dEta, Float_t &dPhi)
TString fInputCaloCellsName
Input cells branch name, if different from default branch.
TGeoHMatrix * fGeomMatrix[22]
Geometry matrices with alignments.
Bool_t AcceptCell(Int_t absID, Bool_t badmap=kTRUE)
Bool_t IsLEDEvent(const Int_t run)
Check if event is LED, is so remove it. Affected LHC11a runs.
AliEMCALRecParam * fRecParam
Reconstruction parameters container.
void SetPositionAlgorithm(Int_t alg)
Int_t fMinEvent
Set a minimum event number, for testing.
TString fGeomName
Name of geometry to use.
Bool_t fOutputAODBranchSet
Set the AOD clusters branch in the input event once.
void SwitchOffL1PhaseInTimeRecalibration()
void PrintParam()
Print clusterization task parameters.
Float_t fTCardCorrMinInduced
Minimum induced energy signal on adjacent cells, sum of induced plus original energy, use same as cell energy clusterization cut.
short Short_t
Definition: External.C:23
TString fOCDBpath
Path with OCDB location.
AliEMCALAfterBurnerUF * fUnfolder
! Unfolding procedure
void SwitchOnDistToBadChannelRecalculation()
Bool_t IsGoodCluster(AliVCluster *cluster, const AliEMCALGeometry *geom, AliVCaloCells *cells, Int_t bc=-1)
Bool_t fGeomMatrixSet
Set geometry matrices only once, for the first event.
Float_t fCentralityBin[2]
Minimum and maximum value of the centrality for the analysis.
Int_t fMaxEvent
Set a maximum event number, for testing.
Bool_t fFillAODCaloCells
Copy calocells to standard branch.
TObjArray * fCaloClusterArr
! CaloClusters array
void RecalculateClusterPosition(const AliEMCALGeometry *geom, AliVCaloCells *cells, AliVCluster *clu)
Double_t fCellTime[fgkNEMCalCells]
Array with cluster time to be passed to digit in case of AODs.
void RecalculateClusterPID(AliVCluster *cluster)
TString fCentralityClass
Name of selected centrality class.
void SetEMCALL1PhaseInTimeRecalibrationForAllSM(const TObjArray *map)
void SwitchOnRejectExoticCluster()
TFile * file
TList with histograms for a given trigger.
void SetExoticCellMinAmplitudeCut(Float_t ma)
Bool_t fAccessOCDB
Need to access info from OCDB (not really)
TObjArray * fClusterArr
! Recpoints array
void FindMatches(AliVEvent *event, TObjArray *clusterArr=0x0, const AliEMCALGeometry *geom=0x0, AliMCEvent *mc=0x0)
unsigned short UShort_t
Definition: External.C:28
void SetExoticCellFractionCut(Float_t f)
Bool_t fRecalibrateWithClusterTime
Use fCellTime to store time of cells in cluster.
void SetEMCALChannelTimeRecalibrationFactors(const TObjArray *map)
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Float_t CorrectClusterEnergyLinearity(AliVCluster *clu)
void FillAODHeader()
Put event header information in standard AOD branch.
Bool_t ClusterContainsBadChannel(const AliEMCALGeometry *geom, const UShort_t *cellList, Int_t nCells)
static const Int_t fgkNEMCalCells
Total number of cells in the calorimeter, 10*48*24 (EMCal) + 4*48*8 (EMCal/DCal 1/3) + 6*32*24 (DCal)...
void CalculateInducedEnergyInTCardCell(Int_t absId, Int_t absIdRef, Int_t sm, Float_t ampRef, Int_t cellCase)
Float_t fTCardCorrInduceEnerFracMin[22]
In case fTCardCorrInduceEnerFracP1 is non null, restrict the minimum fraction of induced energy per S...
void SetEMCALChannelRecalibrationFactors(const TObjArray *map)
Int_t fEMCALNcellsCut
At least an EMCAL cluster with fNCellsCut cells over fEnergyCut.
void SetEMCALChannelStatusMap(const TObjArray *map)
Reclusterize EMCal clusters, put them in a new branch for other following analysis.
Bool_t fTCardCorrCellsNew[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells, that before had no signal.
Bool_t fLoadGeomMatrices
Matrices set from configuration, not get from geometry.root or from ESDs/AODs.
Float_t fEMCALEnergyCut
At least an EMCAL cluster with this energy in the event.
Bool_t fRandomizeTCard
Use random induced energy.
Bool_t fUseAliCentrality
Use the centrality estimator from AliCentrality or AliMultSelection.
Float_t fTCardCorrInduceEnerFracWidth[4][22]
Induced energy loss gauss witdth on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ ...
void RecalibrateCellTime(Int_t absId, Int_t bc, Double_t &time, Bool_t isLGon=kFALSE) const
Int_t fOrgClusterCellId[fgkNEMCalCells]
Array ID of cluster to wich the cell belongs in unmodified clusters.
Float_t fSelectCellMinFrac
Min fraction of cell energy after unfolding cut.
Bool_t fRemapMCLabelForAODs
Remap AOD cells MC label. Needed in old AOD productions.