AliPhysics  857879d (857879d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
44 // --- EMCAL
45 #include "AliEMCALAfterBurnerUF.h"
46 #include "AliEMCALGeometry.h"
47 #include "AliEMCALClusterizerNxN.h"
48 #include "AliEMCALClusterizerv1.h"
49 #include "AliEMCALClusterizerv2.h"
50 #include "AliEMCALRecPoint.h"
51 #include "AliEMCALDigit.h"
52 
54 
58 
59 //______________________________________________________________________________
60 // Constructor.
61 //______________________________________________________________________________
63 : AliAnalysisTaskSE(name)
64 , fEvent(0)
65 , fGeom(0), fGeomName("")
66 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
67 , fOCDBpath(""), fAccessOCDB(kFALSE)
68 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
69 , fRecParam(0), fClusterizer(0)
70 , fUnfolder(0), fJustUnfold(kFALSE)
71 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
72 , fFillAODFile(kFALSE), fFillAODHeader(0)
73 , fFillAODCaloCells(0), fRun(-1)
74 , fRecoUtils(0), fConfigName("")
75 , fOrgClusterCellId()
76 , fCellLabels(), fCellSecondLabels(), fCellTime()
77 , fCellMatchdEta(), fCellMatchdPhi()
78 , fRecalibrateWithClusterTime(0)
79 , fMaxEvent(0), fDoTrackMatching(kFALSE)
80 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
81 , fRejectBelowThreshold(kFALSE)
82 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
83 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
84 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
85 , fConstantTimeShift(0)
86 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
87 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
88 , fSetCellMCLabelFromCluster(0)
89 , fSetCellMCLabelFromEdepFrac(0)
90 , fRemapMCLabelForAODs(0)
91 , fInputFromFilter(0)
92 {
93  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = 0;
94 
95  ResetArrays();
96 
97  fCentralityBin[0] = fCentralityBin[1]=-1;
98 }
99 
100 //______________________________________________________________
102 //______________________________________________________________
104 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
105 , fEvent(0)
106 , fGeom(0), fGeomName("")
107 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
108 , fOCDBpath(""), fAccessOCDB(kFALSE)
109 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
110 , fRecParam(0), fClusterizer(0)
111 , fUnfolder(0), fJustUnfold(kFALSE)
112 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
113 , fFillAODFile(kFALSE), fFillAODHeader(0)
114 , fFillAODCaloCells(0), fRun(-1)
115 , fRecoUtils(0), fConfigName("")
116 , fOrgClusterCellId()
117 , fCellLabels(), fCellSecondLabels(), fCellTime()
118 , fCellMatchdEta(), fCellMatchdPhi()
119 , fRecalibrateWithClusterTime(0)
120 , fMaxEvent(0), fDoTrackMatching(kFALSE)
121 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
122 , fRejectBelowThreshold(kFALSE)
123 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
124 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
125 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
126 , fConstantTimeShift(0)
127 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
128 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
129 , fSetCellMCLabelFromCluster(0)
130 , fSetCellMCLabelFromEdepFrac(0)
131 , fRemapMCLabelForAODs(0)
132 , fInputFromFilter(0)
133 {
134  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = 0;
135 
136  ResetArrays();
137 
138  fCentralityBin[0] = fCentralityBin[1]=-1;
139 }
140 
141 //_______________________________________________________________
143 //_______________________________________________________________
145 {
146  if (fDigitsArr)
147  {
148  fDigitsArr->Clear("C");
149  delete fDigitsArr;
150  }
151 
152  if (fClusterArr)
153  {
154  fClusterArr->Delete();
155  delete fClusterArr;
156  }
157 
158  if (fCaloClusterArr)
159  {
160  fCaloClusterArr->Delete();
161  delete fCaloClusterArr;
162  }
163 
164  if(fClusterizer) delete fClusterizer;
165  if(fUnfolder) delete fUnfolder;
166  if(fRecoUtils) delete fRecoUtils;
167 }
168 
169 //_______________________________________________________
173 //_______________________________________________________
175 {
176  if(!fSelectEMCALEvent) return kTRUE; // accept
177 
178  if(fEMCALEnergyCut <= 0) return kTRUE; // accept
179 
180  Int_t nCluster = fEvent -> GetNumberOfCaloClusters();
181  AliVCaloCells * caloCell = fEvent -> GetEMCALCells();
182  Int_t bc = fEvent -> GetBunchCrossNumber();
183 
184  for(Int_t icalo = 0; icalo < nCluster; icalo++)
185  {
186  AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
187 
188  if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
189  fRecoUtils->IsGoodCluster(clus,fGeom,caloCell,bc))
190  {
191 
192  AliDebug(1, Form("Accept : E %2.2f > %2.2f, nCells %d > %d",
193  clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut));
194 
195  return kTRUE;
196  }
197 
198  }// loop
199 
200  AliDebug(1,"Reject");
201 
202  return kFALSE;
203 }
204 
205 //_______________________________________________
208 //_______________________________________________
210 {
211  //Set it only once
212  if(fOADBSet) return ;
213 
214  Int_t runnumber = InputEvent()->GetRunNumber() ;
215  TString pass = GetPass();
216 
217  AliInfo(Form("Get AODB parameters from EMCAL in %s for run %d, and <%s>",fOADBFilePath.Data(),runnumber,pass.Data()));
218 
219  Int_t nSM = fGeom->GetNumberOfSuperModules();
220 
221  // Bad map
222  if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
223  {
224  AliOADBContainer *contBC=new AliOADBContainer("");
225  contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
226 
227  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
228 
229  if(arrayBC)
230  {
231  fRecoUtils->SwitchOnDistToBadChannelRecalculation();
232  AliInfo("Remove EMCAL bad cells");
233 
234  for (Int_t i=0; i < nSM; ++i)
235  {
236  TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
237 
238  if (hbm)
239  delete hbm;
240 
241  hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
242 
243  if (!hbm)
244  {
245  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
246  continue;
247  }
248 
249  hbm->SetDirectory(0);
250  fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
251 
252  } // loop
253  } else AliInfo("Do NOT remove EMCAL bad channels"); // run array
254 
255  delete contBC;
256  } // Remove bad
257 
258  // Energy Recalibration
259  if(fRecoUtils->IsRecalibrationOn())
260  {
261  AliOADBContainer *contRF=new AliOADBContainer("");
262 
263  contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
264 
265  TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
266 
267  if(recal)
268  {
269  TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
270 
271  if(recalpass)
272  {
273  TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
274 
275  if(recalib)
276  {
277  AliInfo("Recalibrate EMCAL");
278  for (Int_t i=0; i<nSM; ++i)
279  {
280  TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
281 
282  if (h)
283  delete h;
284 
285  h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
286 
287  if (!h)
288  {
289  AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
290  continue;
291  }
292 
293  h->SetDirectory(0);
294 
295  fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
296  } // SM loop
297  } else AliInfo("Do NOT recalibrate EMCAL, no params object array"); // array ok
298  } else AliInfo("Do NOT recalibrate EMCAL, no params for pass"); // array pass ok
299  } else AliInfo("Do NOT recalibrate EMCAL, no params for run"); // run number array ok
300 
301  delete contRF;
302  } // Recalibration on
303 
304  // Energy Recalibration, apply on top of previous calibration factors
305  if(fRecoUtils->IsRunDepRecalibrationOn())
306  {
307  AliOADBContainer *contRFTD=new AliOADBContainer("");
308 
309  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
310 
311  TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
312 
313  //If it did not exist for this run, get closes one
314  if (!htd)
315  {
316  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
317  // let's get the closest runnumber instead then..
318  Int_t lower = 0;
319  Int_t ic = 0;
320  Int_t maxEntry = contRFTD->GetNumberOfEntries();
321 
322  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
323  lower = ic;
324  ic++;
325  }
326 
327  Int_t closest = lower;
328  if ( (ic<maxEntry) &&
329  (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
330  closest = ic;
331  }
332 
333  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
334  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
335  }
336 
337  if(htd)
338  {
339  AliInfo("Recalibrate (Temperature) EMCAL");
340 
341  for (Int_t ism=0; ism<nSM; ++ism)
342  {
343  for (Int_t icol=0; icol<48; ++icol)
344  {
345  for (Int_t irow=0; irow<24; ++irow)
346  {
347  Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
348 
349  Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
350  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
351  //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
352  // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
353  fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
354  } // columns
355  } // rows
356  } // SM loop
357  } else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
358 
359  delete contRFTD;
360  } // Run by Run T calibration
361 
362  // Time Recalibration
363  if(fRecoUtils->IsTimeRecalibrationOn())
364  {
365  AliOADBContainer *contTRF=new AliOADBContainer("");
366 
367  contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
368 
369  TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
370 
371  if(trecal)
372  {
373  // pass number should be pass1 except on Run1 and special cases
374  TString passM = pass;
375  if ( pass=="spc_calo" ) passM = "pass3";
376  if ( runnumber > 209121 ) passM = "pass1"; // run2 periods
377  if ( pass == "muon_calo_pass1" && runnumber > 209121 && runnumber < 244284 )
378  passM = "pass0";//period LHC15a-m
379 
380  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
381 
382  if(trecalpass)
383  {
384  AliInfo("Time Recalibrate EMCAL");
385  for (Int_t ibc = 0; ibc < 4; ++ibc)
386  {
387  TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
388 
389  if (h)
390  delete h;
391 
392  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
393 
394  if (!h)
395  {
396  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
397  continue;
398  }
399 
400  h->SetDirectory(0);
401 
402  fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
403  } // bunch crossing loop
404  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
405  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
406 
407  delete contTRF;
408  } // Time recalibration on
409 
410  // L1 Phase Time Recalibration
411  if( fRecoUtils->IsL1PhaseInTimeRecalibrationOn() )
412  {
413  // init default maps first
414  if (!fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationArray())
415  fRecoUtils->InitEMCALL1PhaseInTimeRecalibration() ;
416 
417  AliOADBContainer *contBC = new AliOADBContainer("");
418 
419  TFile *timeFile=new TFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"read");
420  if (!timeFile || timeFile->IsZombie())
421  {
422  AliFatal(Form("EMCALTimeL1PhaseCalib.root was not found in the path provided: %s",fOADBFilePath.Data()));
423  return ;
424  }
425 
426  if (timeFile) delete timeFile;
427 
428  contBC->InitFromFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"AliEMCALTimeL1PhaseCalib");
429 
430  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
431  if (!arrayBC)
432  {
433  AliError(Form("No external L1 phase in time calibration set for run number: %d", runnumber));
434  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
435  }
436  else
437  {
438  // Only 1 L1 phase correction possible, except special cases
439  TString pass2 = "pass1";
440 
441  if ( pass=="muon_calo_pass1" && runnumber > 209121 && runnumber < 244284 )
442  pass2 = "pass0"; // period LHC15a-m
443 
444  TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(pass2);
445  if (!arrayBCpass)
446  {
447  AliError(Form("No external L1 phase in time calibration set for: %d -%s", runnumber,pass2.Data()));
448  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
449  }
450  else AliInfo("Recalibrate L1 Phase time");
451 
452  if(arrayBCpass)
453  {
454  if ( DebugLevel()>0 ) arrayBCpass->Print();
455 
456  TH1C *h = fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationForAllSM();
457  if (h) delete h;
458 
459  h = (TH1C*)arrayBCpass->FindObject(Form("h%d",runnumber));
460 
461  if (!h)
462  {
463  AliFatal(Form("There is no calibration histogram h%d for this run",runnumber));
464  return;
465  }
466 
467  h->SetDirectory(0);
468  fRecoUtils->SetEMCALL1PhaseInTimeRecalibrationForAllSM(h);
469  }
470  }
471 
472  delete contBC;
473  } // L1 Phase Time Recalibration
474 
475  // Parameters already set once, so do not it again
476  fOADBSet = kTRUE;
477 }
478 
479 //_________________________________________________
482 //_________________________________________________
484 {
485  fEvent = InputEvent();
486  if (!fEvent)
487  {
488  Warning("AccessOCDB","Event not available!!!");
489  return kFALSE;
490  }
491 
492  if (fEvent->GetRunNumber()==fRun)
493  return kTRUE;
494  fRun = fEvent->GetRunNumber();
495 
496  AliDebug(1,"Begin");
497 
498  AliCDBManager *cdb = AliCDBManager::Instance();
499 
500 
501  if (fOCDBpath.Length())
502  {
503  cdb->SetDefaultStorage(fOCDBpath.Data());
504  AliInfo(Form("Default storage %s",fOCDBpath.Data()));
505  }
506 
507  cdb->SetRun(fEvent->GetRunNumber());
508 
509  //
510  // EMCAL from RAW OCDB
511  if (fOCDBpath.Contains("alien:"))
512  {
513  cdb->SetSpecificStorage("EMCAL/Calib/Data","raw://");
514  cdb->SetSpecificStorage("EMCAL/Calib/Time","raw://");
515  cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","raw://");
516  }
517 
518  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
519 
520  return kTRUE;
521 }
522 
523 //_____________________________________________________
528 //_____________________________________________________
530 {
531  fEvent = 0x0;
532 
533  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
534  Int_t eventN = Entry();
535  if(aodIH) eventN = aodIH->GetReadEntry();
536 
537  if (eventN > fMaxEvent)
538  return ;
539 
540  //printf("Clusterizer --- Event %d-- \n",eventN);
541 
542  // Check if input event are embedded events
543  // If so, take output event
544  if (aodIH && aodIH->GetMergeEvents())
545  {
546  fEvent = AODEvent();
547 
548  if(!aodIH->GetMergeEMCALCells())
549  AliFatal("Events merged but not EMCAL cells, check analysis settings!");
550 
551  AliDebug(1,"Use embedded events");
552 
553  AliDebug(1,Form("\t InputEvent N Clusters %d, N Cells %d",
554  InputEvent()->GetNumberOfCaloClusters(),InputEvent()->GetEMCALCells()->GetNumberOfCells()));
555 
556  AliDebug(1,Form("\t MergedEvent N Clusters %d, N Cells %d",
557  aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()));
558 
559 // if(DebugLevel() > 1)
560 // {
561 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
562 // {
563 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
564 // if(sigCluster->IsEMCAL()) AliInfo(Form("\t \t Signal cluster: i %d, E %f",icl,sigCluster->E()));
565 // }
566 // }
567 
568  AliDebug(1,Form("\t OutputEvent N Clusters %d, N Cells %d",
569  AODEvent()->GetNumberOfCaloClusters(), AODEvent()->GetEMCALCells()->GetNumberOfCells()));
570  }
571  else if(fInputFromFilter)
572  {
573  //printf("Get Input From Filtered AOD\n");
574  fEvent = AODEvent();
575  }
576  else
577  {
578  fEvent = InputEvent();
581  }
582 
583  if (!fEvent)
584  {
585  AliError("Event not available");
586  return ;
587  }
588 
589  //Process events if there is a high energy cluster
590  if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
591 
592  //-------------------------------------------------------------------------------------
593  // Reject events if LED was firing, use only for LHC11a data
594  // Reject event if triggered by exotic cell and remove exotic cells if not triggered
595  //-------------------------------------------------------------------------------------
596 
597  if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
598 
599  if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
600 
601  //-------------------------------------------------------------------------------------
602  // Set the cluster array in the event (output or input)
603  //-------------------------------------------------------------------------------------
604 
605  if ( fFillAODFile )
606  {
607  //Magic line to write events to AOD filem put after event rejection
608  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
609  }
610  else if( !fOutputAODBranchSet )
611  {
612  // Create array and put it in the input event, if output AOD not selected, only once
613  InputEvent()->AddObject(fOutputAODBranch);
614  fOutputAODBranchSet = kTRUE;
615  AliInfo(Form("Add AOD branch <%s> to input event",fOutputAODBranchName.Data()));
616  }
617 }
618 
619 //____________________________________________________
624 //____________________________________________________
626 {
627  Int_t nClusters = fEvent->GetNumberOfCaloClusters();
628  Int_t nClustersOrg = 0;
629 
630  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
631  if(aodIH && aodIH->GetEventToMerge()) //Embedding
632  nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
633 
634  ResetArrays();
635 
636  // Loop on original clusters, get MC labels, cluster time (OLD AODs),
637  // or track matching residuals (if matching is not requested)
639  {
640  for (Int_t i = 0; i < nClusters; i++)
641  {
642  AliVCluster *clus = 0;
643  if(aodIH && aodIH->GetEventToMerge()) //Embedding
644  clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
645  else
646  clus = fEvent->GetCaloCluster(i);
647 
648  if(!clus) return;
649 
650  nClustersOrg++;
651 
652  if(!clus->IsEMCAL()) continue;
653 
654  Int_t label = clus->GetLabel();
655  Int_t label2 = -1 ;
656  if (clus->GetNLabels() >=2 ) label2 = clus->GetLabelAt(1) ;
657 
658  //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",
659  // i, clus->GetID(), clus->E(), clus->GetTOF()*1e9, clus->GetNCells(), clus->GetNLabels(), label);
660  //
661  //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++)
662  // printf("%d) Label %d, E dep frac %0.2f; ",
663  // imc, clus->GetLabelAt(imc),clus->GetClusterMCEdepFraction(imc));
664  //if(clus->GetNLabels() > 0) printf("\n");
665 
666  UShort_t * index = clus->GetCellsAbsId() ;
667  for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
668  {
669  fOrgClusterCellId[index[icell]] = i;
670  fCellTime[index[icell]] = clus->GetTOF();
671  fCellMatchdEta[index[icell]] = clus->GetTrackDz();
672  fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
673 
675  {
676  fCellLabels[index[icell]] = label;
677  fCellSecondLabels[index[icell]] = label2;
678  }
679  }
680  }
681  }
682 
683  // Transform CaloCells into Digits
684 
685  Int_t idigit = 0;
686  Int_t id = -1;
687  Float_t amp = -1;
688  Double_t time = -1;
689 
690  AliVCaloCells * cells = fEvent->GetEMCALCells();
691 
692  Int_t bc = InputEvent()->GetBunchCrossNumber();
693 
694  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
695  {
696  // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
697  id = cells->GetCellNumber(icell);
698  Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
699  time-=fConstantTimeShift*1e-9; // only in case of simulations done before 2015
700 
701  // Do not include cells with too low energy, nor exotic cell
702  if( amp < fRecParam->GetMinECut() ||
703  time > fRecParam->GetTimeMax() ||
704  time < fRecParam->GetTimeMin() ) accept = kFALSE;
705 
706  // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
708  {
709  time = fCellTime[id];
710  //printf("cell %d time org %f - ",id, time*1.e9);
711  fRecoUtils->RecalibrateCellTime(id,bc,time);
712  //printf("recal %f\n",time*1.e9);
713  }
714 
715  //Exotic?
716  if (accept && fRecoUtils->IsExoticCell(id,cells,bc))
717  accept = kFALSE;
718 
719  if( !accept )
720  {
721  AliDebug(2,Form("Remove channel absId %d, index %d of %d, amp %f, time %f",
722  id,icell, cells->GetNumberOfCells(), amp, time*1.e9));
723  continue;
724  }
725 
726  //
727  // MC
728  //
729 
730  // Old way
731  Int_t mcLabel = cells->GetMCLabel(icell);
732  Float_t eDep = amp;
733 
734  //printf("--- Selected cell %d--- amp %f\n",id,amp);
735 
736  // New way
737  TArrayI labeArr(0);
738  TArrayF eDepArr(0);
739  Int_t nLabels = 0;
740 
742  {
743  // Old way to recover/set the cell MC label
744  // Only possibility for old Run1 productions
745 
746  if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
747 
748  else if( fSetCellMCLabelFromCluster == 0 &&
750 
751  else mcLabel = -1; // found later
752 
753  // Last parameter of the digit object should be MC deposited energy,
754  // since it is not available in aliroot before year 2016, add just the cell amplitude so that
755  // we give more weight to the MC label of the cell with highest energy in the cluster
756 
757  Float_t efrac = cells->GetEFraction(icell);
758 
759  // When checking the MC of digits, give weight to cells with embedded signal
760  if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
761 
762  eDep *= efrac ;
763  }
764  else // fSetCellMCLabelFromEdepFrac = true
765  {
766  // New way, valid only for MC productions with aliroot > v5-07-21
767  mcLabel = -1;
768 
769  // Map the digit to cell index for later to calculate the cell MC energy deposition map
770  fCellLabels[id] = idigit;
771  //printf("\t absId %d, idigit %d\n",id,idigit);
772 
773  if(fOrgClusterCellId[id] < 0) continue; // index can be negative if noisy cell that did not form cluster
774 
775  AliVCluster *clus = 0;
776  Int_t iclus = fOrgClusterCellId[id];
777 
778  if(iclus < 0)
779  {
780  AliInfo("Negative original cluster index, skip \n");
781  continue;
782  }
783 
784  if(aodIH && aodIH->GetEventToMerge()) //Embedding
785  clus = aodIH->GetEventToMerge()->GetCaloCluster(iclus); //Get clusters directly from embedded signal
786  else
787  clus = fEvent->GetCaloCluster(iclus);
788 
789  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(id, clus, MCEvent(), amp, labeArr, eDepArr);
790  nLabels = labeArr.GetSize();
791 
792  } // cell MC label, new
793 
794  //
795  // Create the digit
796  //
797  if(amp <= 0.01) continue ; // accept if > 10 MeV
798 
799  //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);
800 
801  AliEMCALDigit* digit = new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, eDep);
802 
803  if(nLabels > 0)
804  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
805 
806  idigit++;
807 
808  }
809 
810  fDigitsArr->Sort();
811 
812  //-------------------------------------------------------------------------------------
813  // Do the clusterization
814  //-------------------------------------------------------------------------------------
815 
816  fClusterizer->Digits2Clusters("");
817 
818  //-------------------------------------------------------------------------------------
819  // Transform the recpoints into AliVClusters
820  //-------------------------------------------------------------------------------------
821 
823 
824  if(!fCaloClusterArr)
825  {
826  AliWarning(Form("No array with CaloClusters, input RecPoints entries %d",fClusterArr->GetEntriesFast()));
827  return;
828  }
829 
830  AliDebug(1,Form("N clusters: before recluster %d, after recluster %d, recpoints %d",
831  nClustersOrg, fCaloClusterArr->GetEntriesFast(),fClusterArr->GetEntriesFast()));
832 
833 // if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
834 // {
835 // AliInfo("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
836 // fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
837 // fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
838 // }
839 }
840 
841 
842 //_____________________________________________________
844 //_____________________________________________________
846 {
847  AliVCaloCells *cells = fEvent->GetEMCALCells();
848  Double_t cellAmplitude = 0;
849  Double_t cellTime = 0;
850  Short_t cellNumber = 0;
851  Int_t cellMCLabel = 0;
852  Double_t cellEFrac = 0;
853  Int_t nClustersOrg = 0;
854 
855  // Fill the array with the EMCAL clusters, copy them
856  for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
857  {
858  AliVCluster *clus = fEvent->GetCaloCluster(i);
859  if(clus->IsEMCAL())
860  {
861  //recalibrate/remove bad channels/etc if requested
862  if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
863  {
864  continue;
865  }
866 
867  if(fRecoUtils->IsRecalibrationOn())
868  {
869  //Calibrate cluster
870  fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
871 
872  //CalibrateCells
873  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
874  {
875  if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
876  break;
877 
878  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
879  fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
880  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
881 
882  //Do not include bad channels found in analysis?
883  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
884  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
885  continue;
886 
887  cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
888 
889  }// cells loop
890  }// recalibrate
891 
892  //Cast to ESD or AOD, needed to create the cluster array
893  AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
894  AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
895 
896  if (esdCluster)
897  {
898  fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
899  }//ESD
900  else if(aodCluster)
901  {
902  fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
903  }//AOD
904  else
905  AliWarning("Wrong CaloCluster type?");
906 
907  nClustersOrg++;
908  }
909  }
910 
911  //Do the unfolding
912  fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
913 
914  //CLEAN-UP
915  fUnfolder->Clear();
916 }
917 
918 //_____________________________________________________
920 //_____________________________________________________
922 {
923  AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
924  Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
925 
926  AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
927  aodEMcells.CreateContainer(nEMcell);
928  aodEMcells.SetType(AliVCaloCells::kEMCALCell);
929  Double_t calibFactor = 1.;
930  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
931  {
932  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
933  fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
934  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
935 
936  if(fRecoUtils->IsRecalibrationOn())
937  {
938  calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
939  }
940 
941  if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
942  { //Channel is not declared as bad
943  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
944  eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
945  }
946  else
947  {
948  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
949  }
950  }
951 
952  aodEMcells.Sort();
953 }
954 
955 //__________________________________________________
957 //__________________________________________________
959 {
960  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
961  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
962 
963  Double_t pos[3] ;
964  Double_t covVtx[6];
965  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
966 
967  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
968  if(!header) AliFatal("Not a standard AOD");
969  header->SetRunNumber(fEvent->GetRunNumber());
970 
971  if(esdevent)
972  {
973  TTree* tree = fInputHandler->GetTree();
974  if (tree)
975  {
976  TFile* file = tree->GetCurrentFile();
977  if (file) header->SetESDFileName(file->GetName());
978  }
979  }
980  else if (aodevent) {
981  AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
982  if(!aodheader) AliFatal("Not a standard AOD");
983  header->SetESDFileName(aodheader->GetESDFileName());
984  }
985 
986  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
987  header->SetOrbitNumber(fEvent->GetOrbitNumber());
988  header->SetPeriodNumber(fEvent->GetPeriodNumber());
989  header->SetEventType(fEvent->GetEventType());
990 
991  // Centrality
993  {
994  if(fEvent->GetCentrality())
995  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
996  else
997  header->SetCentrality(0);
998  }
999 
1000  // Trigger
1001  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
1002 
1003  header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
1004 
1005  header->SetTriggerMask(fEvent->GetTriggerMask());
1006  header->SetTriggerCluster(fEvent->GetTriggerCluster());
1007 
1008  if (esdevent)
1009  {
1010  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
1011  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
1012  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
1013  }
1014  else if (aodevent)
1015  {
1016  header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
1017  header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
1018  header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
1019  }
1020 
1021  header->SetMagneticField(fEvent->GetMagneticField());
1022  //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
1023 
1024  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
1025  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
1026  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
1027  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
1028  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
1029 
1030  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
1031  Float_t diamcov[3];
1032  fEvent->GetDiamondCovXY(diamcov);
1033  header->SetDiamond(diamxy,diamcov);
1034  if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
1035  else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
1036 
1037  //
1038  Int_t nVertices = 1 ;/* = prim. vtx*/;
1039  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1040 
1041  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1042 
1043  // Access to the AOD container of vertices
1044  TClonesArray &vertices = *(AODEvent()->GetVertices());
1045  Int_t jVertices=0;
1046 
1047  // Add primary vertex. The primary tracks will be defined
1048  // after the loops on the composite objects (V0, cascades, kinks)
1049  fEvent->GetPrimaryVertex()->GetXYZ(pos);
1050  Float_t chi = 0;
1051  if (esdevent)
1052  {
1053  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1054  chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
1055  }
1056  else if (aodevent)
1057  {
1058  aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1059  chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
1060  }
1061 
1062  AliAODVertex * primary = new(vertices[jVertices++])
1063  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
1064  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
1065  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
1066 }
1067 
1068 //___________________________________________________________
1072 //___________________________________________________________
1074 {
1075  // First recalculate track-matching for the new clusters
1076  if(fDoTrackMatching)
1077  {
1078  fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom,MCEvent());
1079  }
1080 
1081  // Put the new clusters in the AOD list
1082 
1083  Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
1084 
1085  for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
1086  {
1087  AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1088 
1089  newCluster->SetID(i);
1090 
1091  // Correct cluster energy non linearity
1092  newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
1093 
1094  //Add matched track
1095  if(fDoTrackMatching)
1096  {
1097  Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1098  if(trackIndex >= 0)
1099  {
1100  newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1101  AliDebug(2,Form("Matched Track index %d to new cluster %d",trackIndex,i));
1102  }
1103 
1104  Float_t dR = 999., dZ = 999.;
1105  fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dZ,dR);
1106  newCluster->SetTrackDistance(dR,dZ);
1107  }
1108  else
1109  {// Assign previously assigned matched track in reco, very very rough
1110  Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1111  newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1112  }
1113 
1114  //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1115  //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1116  //printf("\n");
1117  //
1118  //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",
1119  // i, newCluster->GetID(), newCluster->E(), newCluster->GetTOF()*1e9, newCluster->GetNCells(), newCluster->GetNLabels(), newCluster->GetLabel());
1120  //
1121  //for(Int_t imc = 0; imc < newCluster->GetNLabels(); imc++)
1122  // printf("%d) Label %d, E dep frac %0.2f; ",
1123  // imc,newCluster->GetLabelAt(imc),newCluster->GetClusterMCEdepFraction(imc));
1124  //if(newCluster->GetNLabels() > 0) printf("\n");
1125 
1126  // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1127  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1128 
1129  new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
1130 
1131  AliDebug(2,Form("New cluster %d of %d, energy %f, mc label %d",
1132  newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()));
1133 
1134  } // cluster loop
1135 
1136  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1137 }
1138 
1139 
1140 //________________________________________________________________
1142 //________________________________________________________________
1144 {
1145  if(fUseAliCentrality)
1146  {
1147  if(GetCentrality()) return GetCentrality()->GetCentralityPercentile(fCentralityClass) ;
1148  else return -1. ;
1149  }
1150  else
1151  {
1152  if(GetMultSelCen()) return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass,kTRUE) ;
1153  else return -1. ;
1154  }
1155 
1156  return -1.;
1157 }
1158 
1159 //_______________________________________________
1161 //_______________________________________________
1163 {
1164  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1165  {
1166  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null");
1167  return TString("");
1168  }
1169 
1170  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1171  {
1172  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null");
1173  return TString("");
1174  }
1175 
1176  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1177  if (pass.Contains("ass1")) return TString("pass1");
1178  else if (pass.Contains("ass2")) return TString("pass2");
1179  else if (pass.Contains("ass3")) return TString("pass3");
1180  else if (pass.Contains("ass4")) return TString("pass4");
1181  else if (pass.Contains("ass5")) return TString("pass5");
1182  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1183  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1184  {
1185  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1186  return TString("pass1");
1187  }
1188  else if (pass.Contains("LHC14a1a"))
1189  {
1190  AliInfo("Enable EMCal energy calibration for this MC production!!");
1191 
1192  fRecoUtils->SwitchOnRecalibration();
1193 
1194  return TString("LHC14a1a");
1195  }
1196 
1197  // No condition fullfilled, give a default value
1198  AliInfo("Pass number string not found");
1199 
1200  return TString("");
1201 }
1202 
1203 //_________________________________________
1206 //_________________________________________
1208 {
1209  if(fDebug >=0) (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1210 
1211  fOADBSet = kFALSE;
1212  if(fOADBFilePath == "") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
1213 
1214  fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1215 
1216  if(!fRecParam) fRecParam = new AliEMCALRecParam;
1217  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1218 
1219  if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1220  if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1221  if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1222  fRejectBelowThreshold = kFALSE;
1223 
1224  // Centrality
1225  if(fCentralityClass == "") fCentralityClass = "V0M";
1226 
1227  if (fOCDBpath == "") fOCDBpath = "raw://" ;
1228  if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1229 
1230  if(gROOT->LoadMacro(fConfigName) >=0)
1231  {
1232  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
1233  AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1234  fGeomName = clus->fGeomName;
1236  fOCDBpath = clus->fOCDBpath;
1237  fAccessOCDB = clus->fAccessOCDB;
1238  fRecParam = clus->fRecParam;
1239  fJustUnfold = clus->fJustUnfold;
1240  fFillAODFile = clus->fFillAODFile;
1241  fRecoUtils = clus->fRecoUtils;
1242  fConfigName = clus->fConfigName;
1243  fMaxEvent = clus->fMaxEvent;
1246  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1248  fCentralityBin[0] = clus->fCentralityBin[0];
1249  fCentralityBin[1] = clus->fCentralityBin[1];
1251  }
1252 }
1253 
1254 //_______________________________________________________
1257 //_______________________________________________________
1259 {
1260  if (fJustUnfold)
1261  {
1262  // init the unfolding afterburner
1263  delete fUnfolder;
1264  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1265  return;
1266  }
1267 
1268  //First init the clusterizer
1269  delete fClusterizer;
1270  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1271  fClusterizer = new AliEMCALClusterizerv1 (fGeom);
1272  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1273  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1274  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1275  {
1276  fClusterizer = new AliEMCALClusterizerNxN(fGeom);
1277  fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1278  fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1279  }
1280  else
1281  {
1282  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1283  }
1284 
1285  //Now set the parameters
1286  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1287  fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1288  fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1289  fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1290  fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1291  fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1292  fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1293  fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1294  fClusterizer->SetTimeCalibration ( fRecParam->IsTimeCalibrationOn() );
1295  fClusterizer->SetInputCalibrated ( kTRUE );
1296  fClusterizer->SetJustClusters ( kTRUE );
1297 
1298  // Initialize the cluster rec points and digits arrays and get them.
1299  fClusterizer->SetOutput(0);
1300  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1301  fDigitsArr = fClusterizer->GetDigits();
1302 
1303  //In case of unfolding after clusterization is requested, set the corresponding parameters
1304  if(fRecParam->GetUnfold())
1305  {
1306  Int_t i=0;
1307  for (i = 0; i < 8; i++)
1308  {
1309  fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1310  }//end of loop over parameters
1311 
1312  for (i = 0; i < 3; i++)
1313  {
1314  fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1315  fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1316  }//end of loop over parameters
1317  fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1318  fClusterizer->InitClusterUnfolding();
1319 
1320  }// to unfold
1321 }
1322 
1323 //_________________________________________________
1327 //_________________________________________________
1329 {
1330  if(fGeomMatrixSet) return;
1331 
1332  Int_t runnumber = InputEvent()->GetRunNumber() ;
1333  if (!fGeom)
1334  {
1335  if(fGeomName=="")
1336  {
1337  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(runnumber);
1338  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fGeom->GetName(),runnumber));
1339  }
1340  else
1341  {
1342  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1343  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeomName.Data()));
1344  }
1345 
1346  // Init geometry, I do not like much to do it like this ...
1347  if(fImportGeometryFromFile && !gGeoManager)
1348  {
1349  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1350  {
1351  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1352  if (runnumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1353  else if (runnumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1354  else if (runnumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1355  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >=2015
1356  }
1357 
1358  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1359 
1360  TGeoManager::Import(fImportGeometryFilePath) ;
1361  }
1362 
1363  AliDebug(1,Form("Init for run=%d",runnumber));
1364  if (!gGeoManager) AliDebug(1,"Careful!, gGeoManager not loaded, load misalign matrices");
1365  } // geometry pointer did not exist before
1366 
1367  if(fLoadGeomMatrices)
1368  {
1369  AliInfo("Load user defined EMCAL geometry matrices");
1370 
1371  // OADB if available
1372  AliOADBContainer emcGeoMat("AliEMCALgeo");
1373  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1374  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1375 
1376  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1377  {
1378 
1379  if (!fGeomMatrix[mod]) // Get it from OADB
1380  {
1381  AliDebug(2,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
1382  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1383 
1384  fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1385  }
1386 
1387  if(fGeomMatrix[mod])
1388  {
1389  if(DebugLevel() > 1)
1390  fGeomMatrix[mod]->Print();
1391 
1392  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1393  }
1394  else if(gGeoManager)
1395  {
1396  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
1397  fGeom->SetMisalMatrix(fGeom->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
1398  }
1399  else
1400  {
1401  AliError(Form("Alignment atrix for SM %d is not available",mod));
1402  }
1403 
1404  fGeomMatrixSet=kTRUE;
1405 
1406  }//SM loop
1407  }//Load matrices
1408  else if(!gGeoManager)
1409  {
1410  AliInfo("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1411  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1412  if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1413  {
1414  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
1415  }//AOD
1416  else
1417  {
1418  AliDebug(1,"Load Misaligned matrices");
1419 
1420  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1421 
1422  if(!esd)
1423  {
1424  AliError("This event does not contain ESDs?");
1425  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1426  return;
1427  }
1428 
1429  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1430  {
1431  if(DebugLevel() > 1)
1432  esd->GetEMCALMatrix(mod)->Print();
1433 
1434  if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1435 
1436  }
1437 
1438  fGeomMatrixSet=kTRUE;
1439 
1440  } // ESD
1441  } // Load matrices from Data
1442 }
1443 
1444 //____________________________________________________
1449 //____________________________________________________
1451 {
1452  if(!fRemoveExoticEvents) return kFALSE;
1453 
1454  // Loop on cells
1455 
1456  AliVCaloCells * cells = fEvent->GetEMCALCells();
1457  Float_t totCellE = 0;
1458  Int_t bc = InputEvent()->GetBunchCrossNumber();
1459 
1460  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1461  {
1462  Float_t ecell = 0 ;
1463  Double_t tcell = 0 ;
1464 
1465  Int_t absID = cells->GetCellNumber(icell);
1466  Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1467  tcell-=fConstantTimeShift*1e-9;// Only for MC simulations done before 2015
1468 
1469  if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1470  }
1471 
1472  // TString triggerclasses = event->GetFiredTriggerClasses();
1473  // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1474  // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1475  // return;
1476  //
1477 
1478  //printf("TotE cell %f\n",totCellE);
1479  if(totCellE < 1) return kTRUE;
1480 
1481  return kFALSE;
1482 }
1483 
1484 //________________________________________________________________
1486 //________________________________________________________________
1488 {
1489  if(!fRemoveLEDEvents) return kFALSE;
1490 
1491  //check events of LHC11a period
1492  if(run < 146858 || run > 146860) return kFALSE ;
1493 
1494  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1495  Int_t ncellsSM3 = 0;
1496  AliVCaloCells * cells = fEvent->GetEMCALCells();
1497  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1498  {
1499  if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1500  }
1501 
1502  TString triggerclasses = fEvent->GetFiredTriggerClasses();
1503 
1504  Int_t ncellcut = 21;
1505  if(triggerclasses.Contains("EMC")) ncellcut = 35;
1506 
1507  if( ncellsSM3 >= ncellcut)
1508  {
1509  AliInfo(Form("Reject event %d with ncells in SM3 %d",(Int_t)Entry(),ncellsSM3));
1510  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1511  return kTRUE;
1512  }
1513 
1514  return kFALSE;
1515 }
1516 
1517 //_______________________________________________________
1521 //_______________________________________________________
1523 {
1524  Int_t j = 0;
1525  for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1526  {
1527  AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1528 
1529  const Int_t ncells = recPoint->GetMultiplicity();
1530  Int_t ncellsTrue = 0;
1531 
1532  if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1533 
1534  // cells and their amplitude fractions
1535  UShort_t absIds[ncells];
1536  Double32_t ratios[ncells];
1537 
1538  //For later check embedding
1539  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1540 
1541  Float_t clusterE = 0;
1542  for (Int_t c = 0; c < ncells; c++)
1543  {
1544  AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1545 
1546  absIds[ncellsTrue] = digit->GetId();
1547  ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1548 
1549  // In case of unfolding, remove digits with unfolded energy too low
1550  if(fSelectCell)
1551  {
1552  if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1553  {
1554 
1555  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
1556  recPoint->GetEnergiesList()[c],digit->GetAmplitude()));
1557  continue;
1558  } // if cuts
1559  }// Select cells
1560 
1561  //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1562  clusterE +=recPoint->GetEnergiesList()[c];
1563 
1564  // In case of embedding, fill ratio with amount of signal,
1565  if (aodIH && aodIH->GetMergeEvents())
1566  {
1567  //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1568  AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1569  AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1570 
1571  Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1572  //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1573  Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1574  //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1575 
1576  if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1577  //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1578 
1579  }//Embedding
1580 
1581  ncellsTrue++;
1582 
1583  }// cluster cell loop
1584 
1585  if (ncellsTrue < 1)
1586  {
1587  AliDebug(2,Form("Skipping cluster with no cells avobe threshold E = %f, ncells %d",
1588  recPoint->GetEnergy(), ncells));
1589  continue;
1590  }
1591 
1592  //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1593 
1594  if(clusterE < fRecParam->GetClusteringThreshold())
1595  {
1596  AliDebug(2,Form("Remove cluster with energy below seed threshold %f",clusterE));
1597  continue;
1598  }
1599 
1600  TVector3 gpos;
1601  Float_t g[3];
1602 
1603  // calculate new cluster position
1604 
1605  recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1606  recPoint->GetGlobalPosition(gpos);
1607  gpos.GetXYZ(g);
1608 
1609  // create a new cluster
1610 
1611  (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1612  AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1613  j++;
1614  clus->SetType(AliVCluster::kEMCALClusterv1);
1615  clus->SetE(clusterE);
1616  clus->SetPosition(g);
1617  clus->SetNCells(ncellsTrue);
1618  clus->SetCellsAbsId(absIds);
1619  clus->SetCellsAmplitudeFraction(ratios);
1620  clus->SetChi2(-1); //not yet implemented
1621  clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1622  clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1623  clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1624 
1625  if(ncells == ncellsTrue)
1626  {
1627  Float_t elipAxis[2];
1628  recPoint->GetElipsAxis(elipAxis);
1629  clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1630  clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1631  clus->SetDispersion(recPoint->GetDispersion());
1632  }
1633  else if(fSelectCell)
1634  {
1635  // In case some cells rejected, in unfolding case, recalculate
1636  // shower shape parameters and position
1637  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
1638 
1639  AliVCaloCells* cells = 0x0;
1640  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1641  else cells = InputEvent()->GetEMCALCells();
1642 
1643  fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1644  fRecoUtils->RecalculateClusterPID(clus);
1645  fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1646  }
1647 
1648  // MC
1649 
1652  else
1653  {
1654  //
1655  // Normal case, trust what the clusterizer has found
1656  //
1657  Int_t parentMult = 0;
1658  Int_t *parentList = recPoint->GetParents(parentMult);
1659  Float_t *parentListDE = recPoint->GetParentsDE(); // deposited energy
1660 
1661  clus->SetLabel(parentList, parentMult);
1663  clus->SetClusterMCEdepFractionFromEdepArray(parentListDE);
1664 
1665  //
1666  // Set the cell energy deposition fraction map:
1667  //
1668  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
1669  {
1670  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
1671 
1672  // Get the digit that originated this cell cluster
1673  AliVCaloCells* cells = 0x0;
1674  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1675  else cells = InputEvent()->GetEMCALCells();
1676 
1677  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
1678  {
1679  Int_t idigit = fCellLabels[absIds[icell]];
1680 
1681  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
1682 
1683  // Find the 4 MC labels that contributed to the cluster and their
1684  // deposited energy in the current digit
1685 
1686  mcEdepFracPerCell[icell] = 0; // init
1687 
1688  Int_t nparents = dig->GetNiparent();
1689  if ( nparents > 0 )
1690  {
1691  Int_t digLabel =-1 ;
1692  Float_t edep = 0 ;
1693  Float_t edepTot = 0 ;
1694  Float_t mcEDepFrac[4] = {0,0,0,0};
1695 
1696  // all parents in digit
1697  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
1698  {
1699  digLabel = dig->GetIparent (jndex+1);
1700  edep = dig->GetDEParent(jndex+1);
1701  edepTot += edep;
1702 
1703  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
1704  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
1705  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
1706  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
1707  } // all prarents in digit
1708 
1709  // Divide energy deposit by total deposited energy
1710  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
1711  if(edepTot > 0.01)
1712  {
1713  mcEdepFracPerCell[icell] = clus->PackMCEdepFraction(mcEDepFrac);
1714  }
1715  } // at least one parent label in digit
1716  } // cell in cluster loop
1717 
1718  clus->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
1719 
1720  delete [] mcEdepFracPerCell;
1721 
1722  } // at least one parent in cluster, do the cell primary packing
1723  }
1724  } // recPoints loop
1725 }
1726 
1727 //_____________________________________________________________________
1730 //_____________________________________________________________________
1732 {
1733  if(label < 0) return ;
1734 
1735  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
1736  if(!evt) return ;
1737 
1738  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1739  if(!arr) return ;
1740 
1741  if(label < arr->GetEntriesFast())
1742  {
1743  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1744  if(!particle) return ;
1745 
1746  if(label == particle->Label()) return ; // label already OK
1747  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
1748  }
1749  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
1750 
1751  // loop on the particles list and check if there is one with the same label
1752  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
1753  {
1754  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1755  if(!particle) continue ;
1756 
1757  if(label == particle->Label())
1758  {
1759  label = ind;
1760  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
1761  return;
1762  }
1763  }
1764 
1765  label = -1;
1766 
1767  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
1768 }
1769 
1770 //________________________________________________
1772 //________________________________________________
1774 {
1775  for(Int_t j = 0; j < fgkNEMCalCells; j++)
1776  {
1777  fOrgClusterCellId[j] = -1;
1778  fCellLabels[j] = -1;
1779  fCellSecondLabels[j] = -1;
1780  fCellTime[j] = 0.;
1781  fCellMatchdEta[j] = -999;
1782  fCellMatchdPhi[j] = -999;
1783  }
1784 }
1785 
1786 //_____________________________________________________________________________________________________
1790 //_____________________________________________________________________________________________________
1792  AliAODCaloCluster * clus)
1793 {
1794  Int_t parentMult = 0;
1795  Int_t *parentList = recPoint->GetParents(parentMult);
1796  clus->SetLabel(parentList, parentMult);
1797 
1798  //Write the second major contributor to each MC cluster.
1799  Int_t iNewLabel ;
1800  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1801  {
1802 
1803  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1804  if(idCell>=0)
1805  {
1806  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
1807  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1808  {
1809  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
1810  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
1811  }
1812  if (iNewLabel == 1)
1813  {
1814  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1815  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1816  {
1817  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1818  }
1819 
1820  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1821  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1822  delete [] newLabelArray;
1823  }
1824  } // positive cell number
1825  }
1826 }
1827 
1828 //___________________________________________________________________________________________________
1832 //___________________________________________________________________________________________________
1834 {
1835  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1836  clArray.Reset();
1837  Int_t nClu = 0;
1838  Int_t nLabTotOrg = 0;
1839  Float_t emax = -1;
1840  Int_t idMax = -1;
1841 
1842  AliVEvent * event = fEvent;
1843 
1844  // In case of embedding the MC signal is in the added event
1845  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1846  if(aodIH && aodIH->GetEventToMerge()) //Embedding
1847  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
1848 
1849 
1850  // Find the clusters that originally had the cells
1851  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1852  {
1853  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1854 
1855  if(idCell>=0)
1856  {
1857  Int_t idCluster = fOrgClusterCellId[idCell];
1858 
1859  Bool_t set = kTRUE;
1860  for(Int_t icl =0; icl < nClu; icl++)
1861  {
1862  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
1863  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
1864  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
1865  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
1866  }
1867  if( set && idCluster >= 0)
1868  {
1869  clArray.SetAt(idCluster,nClu++);
1870  //printf("******** idCluster %d \n",idCluster);
1871  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
1872 
1873  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
1874 
1875  //Search highest E cluster
1876  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1877  //printf("\t E %f\n",clOrg->E());
1878  if(emax < clOrg->E())
1879  {
1880  emax = clOrg->E();
1881  idMax = idCluster;
1882  }
1883  }
1884  }
1885  }// cell loop
1886 
1887  if(nClu==0 || nLabTotOrg == 0)
1888  {
1889  //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());
1890  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
1891  //printf("\n");
1892  }
1893 
1894  // Put the first in the list the cluster with highest energy
1895  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1896  {
1897  Int_t maxIndex = -1;
1898  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1899  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1900  {
1901  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1902  }
1903 
1904  if(firstCluster >=0 && idMax >=0)
1905  {
1906  clArray.SetAt(idMax,0);
1907  clArray.SetAt(firstCluster,maxIndex);
1908  }
1909  }
1910 
1911  // Get the labels list in the original clusters, assign all to the new cluster
1912  TArrayI clMCArray(nLabTotOrg) ;
1913  clMCArray.Reset();
1914 
1915  Int_t nLabTot = 0;
1916  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1917  {
1918  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
1919  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
1920  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1921  Int_t nLab = clOrg->GetNLabels();
1922 
1923  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
1924  {
1925  Int_t lab = clOrg->GetLabelAt(iLab) ;
1926  if(lab>=0)
1927  {
1928  Bool_t set = kTRUE;
1929  //printf("\t \t Set Label %d \n", lab);
1930  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1931  {
1932  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
1933  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
1934  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
1935  }
1936  if( set ) clMCArray.SetAt(lab,nLabTot++);
1937  }
1938  }
1939  }// cluster loop
1940 
1941  // Set the final list of labels
1942 
1943  clus->SetLabel(clMCArray.GetArray(), nLabTot);
1944 
1945 // printf("Final list of labels for new cluster : \n");
1946 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
1947 // {
1948 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
1949 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
1950 // printf(" org %d ",label);
1951 // RemapMCLabelForAODs(label);
1952 // printf(" new %d \n",label);
1953 // }
1954 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
1955 }
1956 
1957 //____________________________________________________________
1958 // Init geometry, create list of output clusters.
1959 //____________________________________________________________
1961 {
1962  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1963 
1964  if(fOutputAODBranchName.Length()==0)
1965  {
1966  fOutputAODBranchName = "newEMCALClustersArray";
1967  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
1968  }
1969 
1971 
1972  if( fFillAODFile )
1973  {
1974  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1975 
1976  //fOutputAODBranch->SetOwner(kFALSE);
1977 
1978  AddAODBranch("TClonesArray", &fOutputAODBranch);
1979  }
1980 }
1981 
1982 //_______________________________________________________
1993 //_______________________________________________________
1995 {
1996  //-------
1997  // Step 1
1998 
1999  // Remove the contents of AOD branch output list set in the previous event
2000  fOutputAODBranch->Clear("C");
2001 
2002  LoadBranches();
2003 
2004  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
2005  // If we need a centrality bin, we select only those events in the corresponding bin.
2006  if( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
2007  {
2008  Float_t cen = GetEventCentrality();
2009  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
2010  }
2011 
2012  // Intermediate array with new clusters : init the array only once or clear from previous event
2013  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
2014  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
2015 
2016  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
2017 
2018  // Get the event, do some checks and settings
2019  CheckAndGetEvent() ;
2020 
2021  if (!fEvent)
2022  {
2023  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
2024  return ;
2025  }
2026 
2027  // Init pointers, geometry, clusterizer, ocdb, aodb
2028 
2029  if(fAccessOCDB) AccessOCDB();
2030  if(fAccessOADB) AccessOADB(); // only once
2031 
2033 
2034  //-------
2035  // Step 2
2036 
2037  // Make clusters
2039  else ClusterizeCells() ;
2040 
2041  //-------
2042  // Step 3
2043 
2045 }
2046 
2047 
void FillAODCaloCells()
Put calo cells in standard branch.
double Double_t
Definition: External.C:58
virtual void ResetArrays()
Reset arrays containing information for all possible cells.
Definition: External.C:236
TString GetPass()
Get or guess pass number/string from path of filename.
Bool_t fFillAODHeader
Copy header to standard branch.
TClonesArray * fDigitsArr
! Digits array
AliEMCALClusterizer * fClusterizer
! EMCAL clusterizer
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.
TCanvas * c
Definition: TestFitELoss.C:172
Int_t fCellSecondLabels[fgkNEMCalCells]
Array with Second MC label to be passed to digit.
void SetClustersMCLabelFrom2SelectedLabels(AliEMCALRecPoint *recPoint, AliAODCaloCluster *clus)
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.
TString fImportGeometryFilePath
path fo geometry.root file
Bool_t fSelectEMCALEvent
Process the event if there is some high energy cluster.
Float_t GetEventCentrality() const
Get centrality/multiplicity percentile.
Bool_t fInputFromFilter
Get the input from AODs from the filter.
TString fOutputAODBranchName
New of output AOD branch.
AliEMCALGeometry * fGeom
EMCAL geometry.
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.
unsigned int UInt_t
Definition: External.C:33
AliEMCALRecoUtils * fRecoUtils
Access to factorized reconstruction algorithms.
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)
Float_t fSelectCellMinE
Min energy cell threshold, after unfolding.
Bool_t fOADBSet
AODB parameters already set.
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 recalulation to speed up analysis in PbPb.
Float_t fCellMatchdPhi[fgkNEMCalCells]
Array with cluster-track dEta.
TGeoHMatrix * fGeomMatrix[22]
Geometry matrices with alignments.
Bool_t IsLEDEvent(const Int_t run)
Check if event is LED, is so remove it. Affected LHC11a runs.
AliEMCALRecParam * fRecParam
Reconstruction parameters container.
TString fGeomName
Name of geometry to use.
Bool_t fOutputAODBranchSet
Set the AOD clusters branch in the input event once.
short Short_t
Definition: External.C:23
TString fOCDBpath
Path with OCDB location.
AliEMCALAfterBurnerUF * fUnfolder
! Unfolding procedure
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.
Bool_t fFillAODCaloCells
Copy calocells to standard branch.
TObjArray * fCaloClusterArr
! CaloClusters array
Double_t fCellTime[fgkNEMCalCells]
Array with cluster time to be passed to digit in case of AODs.
TString fCentralityClass
Name of selected centrality class.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
TFile * file
TList with histograms for a given trigger.
Bool_t fAccessOCDB
Need to access info from OCDB (not really)
TObjArray * fClusterArr
! Recpoints array
unsigned short UShort_t
Definition: External.C:28
Bool_t fRecalibrateWithClusterTime
Use fCellTime to store time of cells in cluster.
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void FillAODHeader()
Put event header information in standard AOD branch.
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)...
Int_t fEMCALNcellsCut
At least an EMCAL cluster with fNCellsCut cells over fEnergyCut.
Reclusterize EMCal clusters, put them in a new branch for other following analysis.
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 fUseAliCentrality
Use the centrality estimator from AliCentrality or AliMultSelection.
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.