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