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