AliPhysics  32b88a8 (32b88a8)
 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  TString passM = pass;
374  if(pass=="spc_calo") passM = "pass1";
375  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
376 
377  if(trecalpass)
378  {
379  AliInfo("Time Recalibrate EMCAL");
380  for (Int_t ibc = 0; ibc < 4; ++ibc)
381  {
382  TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
383 
384  if (h)
385  delete h;
386 
387  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
388 
389  if (!h)
390  {
391  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
392  continue;
393  }
394 
395  h->SetDirectory(0);
396 
397  fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
398  } // bunch crossing loop
399  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
400  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
401 
402  delete contTRF;
403  } // Time recalibration on
404 
405  // L1 Phase Time Recalibration
406  if( fRecoUtils->IsL1PhaseInTimeRecalibrationOn() )
407  {
408  // init default maps first
409  if (!fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationArray())
410  fRecoUtils->InitEMCALL1PhaseInTimeRecalibration() ;
411 
412  AliOADBContainer *contBC = new AliOADBContainer("");
413 
414  TFile *timeFile=new TFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"read");
415  if (!timeFile || timeFile->IsZombie())
416  {
417  AliFatal(Form("EMCALTimeL1PhaseCalib.root was not found in the path provided: %s",fOADBFilePath.Data()));
418  return ;
419  }
420 
421  if (timeFile) delete timeFile;
422 
423  contBC->InitFromFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"AliEMCALTimeL1PhaseCalib");
424 
425  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
426  if (!arrayBC)
427  {
428  AliError(Form("No external L1 phase in time calibration set for run number: %d", runnumber));
429  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
430  }
431  else
432  {
433  // Here, it looks for a specific pass
434  TString pass2 = pass;
435  if (pass=="calo_spc") pass2 ="pass1";
436  if (pass=="muon_calo_pass1") pass2 ="pass0";
437  if (pass=="muon_calo_pass2" || pass=="pass2" || pass=="pass3" || pass=="pass4") pass2 ="pass1";
438 
439  TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(pass2);
440  if (!arrayBCpass)
441  {
442  AliError(Form("No external L1 phase in time calibration set for: %d -%s", runnumber,pass2.Data()));
443  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
444  }
445  else AliInfo("Recalibrate L1 Phase time");
446 
447  if(arrayBCpass)
448  {
449  if ( DebugLevel()>0 ) arrayBCpass->Print();
450 
451  TH1C *h = fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationForAllSM();
452  if (h) delete h;
453 
454  h = (TH1C*)arrayBCpass->FindObject(Form("h%d",runnumber));
455 
456  if (!h)
457  {
458  AliFatal(Form("There is no calibration histogram h%d for this run",runnumber));
459  return;
460  }
461 
462  h->SetDirectory(0);
463  fRecoUtils->SetEMCALL1PhaseInTimeRecalibrationForAllSM(h);
464  }
465  }
466 
467  delete contBC;
468  } // L1 Phase Time Recalibration
469 
470  // Parameters already set once, so do not it again
471  fOADBSet = kTRUE;
472 }
473 
474 //_________________________________________________
477 //_________________________________________________
479 {
480  fEvent = InputEvent();
481  if (!fEvent)
482  {
483  Warning("AccessOCDB","Event not available!!!");
484  return kFALSE;
485  }
486 
487  if (fEvent->GetRunNumber()==fRun)
488  return kTRUE;
489  fRun = fEvent->GetRunNumber();
490 
491  AliDebug(1,"Begin");
492 
493  AliCDBManager *cdb = AliCDBManager::Instance();
494 
495 
496  if (fOCDBpath.Length())
497  {
498  cdb->SetDefaultStorage(fOCDBpath.Data());
499  AliInfo(Form("Default storage %s",fOCDBpath.Data()));
500  }
501 
502  cdb->SetRun(fEvent->GetRunNumber());
503 
504  //
505  // EMCAL from RAW OCDB
506  if (fOCDBpath.Contains("alien:"))
507  {
508  cdb->SetSpecificStorage("EMCAL/Calib/Data","raw://");
509  cdb->SetSpecificStorage("EMCAL/Calib/Time","raw://");
510  cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","raw://");
511  }
512 
513  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
514 
515  return kTRUE;
516 }
517 
518 //_____________________________________________________
523 //_____________________________________________________
525 {
526  fEvent = 0x0;
527 
528  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
529  Int_t eventN = Entry();
530  if(aodIH) eventN = aodIH->GetReadEntry();
531 
532  if (eventN > fMaxEvent)
533  return ;
534 
535  //printf("Clusterizer --- Event %d-- \n",eventN);
536 
537  // Check if input event are embedded events
538  // If so, take output event
539  if (aodIH && aodIH->GetMergeEvents())
540  {
541  fEvent = AODEvent();
542 
543  if(!aodIH->GetMergeEMCALCells())
544  AliFatal("Events merged but not EMCAL cells, check analysis settings!");
545 
546  AliDebug(1,"Use embedded events");
547 
548  AliDebug(1,Form("\t InputEvent N Clusters %d, N Cells %d",
549  InputEvent()->GetNumberOfCaloClusters(),InputEvent()->GetEMCALCells()->GetNumberOfCells()));
550 
551  AliDebug(1,Form("\t MergedEvent N Clusters %d, N Cells %d",
552  aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()));
553 
554 // if(DebugLevel() > 1)
555 // {
556 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
557 // {
558 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
559 // if(sigCluster->IsEMCAL()) AliInfo(Form("\t \t Signal cluster: i %d, E %f",icl,sigCluster->E()));
560 // }
561 // }
562 
563  AliDebug(1,Form("\t OutputEvent N Clusters %d, N Cells %d",
564  AODEvent()->GetNumberOfCaloClusters(), AODEvent()->GetEMCALCells()->GetNumberOfCells()));
565  }
566  else if(fInputFromFilter)
567  {
568  //printf("Get Input From Filtered AOD\n");
569  fEvent = AODEvent();
570  }
571  else
572  {
573  fEvent = InputEvent();
576  }
577 
578  if (!fEvent)
579  {
580  AliError("Event not available");
581  return ;
582  }
583 
584  //Process events if there is a high energy cluster
585  if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
586 
587  //-------------------------------------------------------------------------------------
588  // Reject events if LED was firing, use only for LHC11a data
589  // Reject event if triggered by exotic cell and remove exotic cells if not triggered
590  //-------------------------------------------------------------------------------------
591 
592  if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
593 
594  if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
595 
596  //-------------------------------------------------------------------------------------
597  // Set the cluster array in the event (output or input)
598  //-------------------------------------------------------------------------------------
599 
600  if ( fFillAODFile )
601  {
602  //Magic line to write events to AOD filem put after event rejection
603  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
604  }
605  else if( !fOutputAODBranchSet )
606  {
607  // Create array and put it in the input event, if output AOD not selected, only once
608  InputEvent()->AddObject(fOutputAODBranch);
609  fOutputAODBranchSet = kTRUE;
610  AliInfo(Form("Add AOD branch <%s> to input event",fOutputAODBranchName.Data()));
611  }
612 }
613 
614 //____________________________________________________
619 //____________________________________________________
621 {
622  Int_t nClusters = fEvent->GetNumberOfCaloClusters();
623  Int_t nClustersOrg = 0;
624 
625  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
626  if(aodIH && aodIH->GetEventToMerge()) //Embedding
627  nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
628 
629  ResetArrays();
630 
631  // Loop on original clusters, get MC labels, cluster time (OLD AODs),
632  // or track matching residuals (if matching is not requested)
634  {
635  for (Int_t i = 0; i < nClusters; i++)
636  {
637  AliVCluster *clus = 0;
638  if(aodIH && aodIH->GetEventToMerge()) //Embedding
639  clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
640  else
641  clus = fEvent->GetCaloCluster(i);
642 
643  if(!clus) return;
644 
645  if(clus->IsEMCAL())
646  {
647  Int_t label = clus->GetLabel();
648  Int_t label2 = -1 ;
649  if (clus->GetNLabels() >=2 ) label2 = clus->GetLabelAt(1) ;
650 
651  //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",
652  // i, clus->GetID(), clus->E(), clus->GetTOF()*1e9, clus->GetNCells(), clus->GetNLabels(), label);
653  //
654  //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++)
655  // printf("%d) Label %d, E dep frac %0.2f; ",
656  // imc, clus->GetLabelAt(imc),clus->GetClusterMCEdepFraction(imc));
657  //if(clus->GetNLabels() > 0) printf("\n");
658 
659  UShort_t * index = clus->GetCellsAbsId() ;
660  for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
661  {
662  fOrgClusterCellId[index[icell]] = i;
663  fCellTime[index[icell]] = clus->GetTOF();
664  fCellMatchdEta[index[icell]] = clus->GetTrackDz();
665  fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
666 
668  {
669  fCellLabels[index[icell]] = label;
670  fCellSecondLabels[index[icell]] = label2;
671  }
672 
673  Float_t eDepFrac[4];
674  clus->GetCellMCEdepFractionArray(icell,eDepFrac);
675  }
676 
677  nClustersOrg++;
678  }
679  }
680  }
681  // Transform CaloCells into Digits
682 
683  Int_t idigit = 0;
684  Int_t id = -1;
685  Float_t amp = -1;
686  Double_t time = -1;
687 
688  AliVCaloCells * cells = fEvent->GetEMCALCells();
689 
690  Int_t bc = InputEvent()->GetBunchCrossNumber();
691 
692  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
693  {
694  // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
695  id = cells->GetCellNumber(icell);
696  Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
697  time-=fConstantTimeShift*1e-9; // only in case of simulations done before 2015
698 
699  // Do not include cells with too low energy, nor exotic cell
700  if( amp < fRecParam->GetMinECut() ||
701  time > fRecParam->GetTimeMax() ||
702  time < fRecParam->GetTimeMin() ) accept = kFALSE;
703 
704  // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
706  {
707  time = fCellTime[id];
708  //printf("cell %d time org %f - ",id, time*1.e9);
709  fRecoUtils->RecalibrateCellTime(id,bc,time);
710  //printf("recal %f\n",time*1.e9);
711  }
712 
713  //Exotic?
714  if (accept && fRecoUtils->IsExoticCell(id,cells,bc))
715  accept = kFALSE;
716 
717  if( !accept )
718  {
719  AliDebug(2,Form("Remove channel absId %d, index %d of %d, amp %f, time %f",
720  id,icell, cells->GetNumberOfCells(), amp, time*1.e9));
721  continue;
722  }
723 
724  Int_t mcLabel = cells->GetMCLabel(icell);
725 
726  // Old way to recover/set the cell MC label
728  {
729  if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
730 
731  else if( fSetCellMCLabelFromCluster == 0 &&
733 
734  else mcLabel = -1; // found later
735  }
736 
737  // Create the digit, put a fake primary deposited energy to trick the clusterizer
738  // when checking the most likely primary
739 
740  Float_t efrac = cells->GetEFraction(icell);
741 
742  //When checking the MC of digits, give weight to cells with embedded signal
743  if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
744 
745  AliEMCALDigit* digit = new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, amp*efrac);
746 
747  // Last parameter should be MC deposited energy, since it is not available in aliroot before year 2016, add just the cell amplitude so that
748  // we give more weight to the MC label of the cell with highest energy in the cluster
749 
750  // New way, valid only for MC productions with aliroot > v5-07-21
752  {
753  // Map the digit to cell index for later to calculate the cell MC energy deposition map
754  fCellLabels[id] = idigit;
755  //printf("\t absId %d, idigit %d\n",id,idigit);
756 
757  if(fOrgClusterCellId[id] >= 0) // index can be negative if noisy cell that did not form cluster
758  {
759  AliVCluster *clus = 0;
760  Int_t iclus = fOrgClusterCellId[id];
761 
762  if(iclus < 0)
763  {
764  AliInfo("Negative original cluster index, skip \n");
765  continue;
766  }
767 
768  if(aodIH && aodIH->GetEventToMerge()) //Embedding
769  clus = aodIH->GetEventToMerge()->GetCaloCluster(iclus); //Get clusters directly from embedded signal
770  else
771  clus = fEvent->GetCaloCluster(iclus);
772 
773  for(Int_t icluscell=0; icluscell < clus->GetNCells(); icluscell++ )
774  {
775  if(id != clus->GetCellAbsId(icluscell)) continue ;
776 
777  // Get the energy deposition fraction.
778  Float_t eDepFrac[4];
779  clus->GetCellMCEdepFractionArray(icluscell,eDepFrac);
780 
781  // Select the MC label contributing, only if enough energy deposition
782  TArrayI labeArr(0);
783  TArrayF eDepArr(0);
784  Int_t nLabels = 0;
785  for(Int_t imc = 0; imc < 4; imc++)
786  {
787  if(eDepFrac[imc] > 0 && clus->GetNLabels() > imc)
788  {
789  nLabels++;
790 
791  labeArr.Set(nLabels);
792  labeArr.AddAt(clus->GetLabelAt(imc), nLabels-1);
793 
794  eDepArr.Set(nLabels);
795  eDepArr.AddAt(eDepFrac[imc]*amp , nLabels-1);
796  // use as deposited energy a fraction of the simulated energy (smeared and with noise)
797  }
798  }
799 
800  if(nLabels > 0)
801  {
802  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
803  }
804  }
805  }
806  } // cell MC label, new
807 
808  idigit++;
809  }
810 
811  fDigitsArr->Sort();
812 
813  //-------------------------------------------------------------------------------------
814  // Do the clusterization
815  //-------------------------------------------------------------------------------------
816 
817  fClusterizer->Digits2Clusters("");
818 
819  //-------------------------------------------------------------------------------------
820  // Transform the recpoints into AliVClusters
821  //-------------------------------------------------------------------------------------
822 
824 
825  if(!fCaloClusterArr)
826  {
827  AliWarning(Form("No array with CaloClusters, input RecPoints entries %d",fClusterArr->GetEntriesFast()));
828  return;
829  }
830 
831  AliDebug(1,Form("N clusters: before recluster %d, after recluster %d, recpoints %d",
832  nClustersOrg, fCaloClusterArr->GetEntriesFast(),fClusterArr->GetEntriesFast()));
833 
834 // if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
835 // {
836 // AliInfo("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
837 // fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
838 // fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
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);
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);
1662  clus->SetClusterMCEdepFractionFromEdepArray(parentListDE);
1663 
1664  //
1665  // Set the cell energy deposition fraction map:
1666  //
1667  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
1668  {
1669  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
1670 
1671  // Get the digit that originated this cell cluster
1672  AliVCaloCells* cells = 0x0;
1673  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1674  else cells = InputEvent()->GetEMCALCells();
1675 
1676  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
1677  {
1678  Int_t idigit = fCellLabels[absIds[icell]];
1679 
1680  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
1681 
1682  // Find the 4 MC labels that contributed to the cluster and their
1683  // deposited energy in the current digit
1684 
1685  mcEdepFracPerCell[icell] = 0; // init
1686 
1687  Int_t nparents = dig->GetNiparent();
1688  if ( nparents > 0 )
1689  {
1690  Int_t digLabel =-1 ;
1691  Float_t edep = 0 ;
1692  Float_t edepTot = 0 ;
1693  Float_t mcEDepFrac[4] = {0,0,0,0};
1694 
1695  // all parents in digit
1696  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
1697  {
1698  digLabel = dig->GetIparent (jndex+1);
1699  edep = dig->GetDEParent(jndex+1);
1700  edepTot += edep;
1701 
1702  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
1703  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
1704  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
1705  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
1706  } // all prarents in digit
1707 
1708  // Divide energy deposit by total deposited energy
1709  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
1710  if(edepTot > 0.01)
1711  {
1712  mcEdepFracPerCell[icell] = clus->PackMCEdepFraction(mcEDepFrac);
1713  }
1714  } // at least one parent label in digit
1715  } // cell in cluster loop
1716 
1717  clus->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
1718 
1719  delete [] mcEdepFracPerCell;
1720 
1721  } // at least one parent in cluster, do the cell primary packing
1722  }
1723  } // recPoints loop
1724 }
1725 
1726 //_____________________________________________________________________
1729 //_____________________________________________________________________
1731 {
1732  if(label < 0) return ;
1733 
1734  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
1735  if(!evt) return ;
1736 
1737  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1738  if(!arr) return ;
1739 
1740  if(label < arr->GetEntriesFast())
1741  {
1742  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1743  if(!particle) return ;
1744 
1745  if(label == particle->Label()) return ; // label already OK
1746  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
1747  }
1748  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
1749 
1750  // loop on the particles list and check if there is one with the same label
1751  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
1752  {
1753  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1754  if(!particle) continue ;
1755 
1756  if(label == particle->Label())
1757  {
1758  label = ind;
1759  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
1760  return;
1761  }
1762  }
1763 
1764  label = -1;
1765 
1766  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
1767 }
1768 
1769 //________________________________________________
1771 //________________________________________________
1773 {
1774  for(Int_t j = 0; j < fgkNEMCalCells; j++)
1775  {
1776  fOrgClusterCellId[j] = -1;
1777  fCellLabels[j] = -1;
1778  fCellSecondLabels[j] = -1;
1779  fCellTime[j] = 0.;
1780  fCellMatchdEta[j] = -999;
1781  fCellMatchdPhi[j] = -999;
1782  }
1783 }
1784 
1785 //_____________________________________________________________________________________________________
1789 //_____________________________________________________________________________________________________
1791  AliAODCaloCluster * clus)
1792 {
1793  Int_t parentMult = 0;
1794  Int_t *parentList = recPoint->GetParents(parentMult);
1795  clus->SetLabel(parentList, parentMult);
1796 
1797  //Write the second major contributor to each MC cluster.
1798  Int_t iNewLabel ;
1799  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1800  {
1801 
1802  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1803  if(idCell>=0)
1804  {
1805  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
1806  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1807  {
1808  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
1809  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
1810  }
1811  if (iNewLabel == 1)
1812  {
1813  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1814  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1815  {
1816  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1817  }
1818 
1819  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1820  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1821  delete [] newLabelArray;
1822  }
1823  } // positive cell number
1824  }
1825 }
1826 
1827 //___________________________________________________________________________________________________
1831 //___________________________________________________________________________________________________
1833 {
1834  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1835  clArray.Reset();
1836  Int_t nClu = 0;
1837  Int_t nLabTotOrg = 0;
1838  Float_t emax = -1;
1839  Int_t idMax = -1;
1840 
1841  AliVEvent * event = fEvent;
1842 
1843  // In case of embedding the MC signal is in the added event
1844  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1845  if(aodIH && aodIH->GetEventToMerge()) //Embedding
1846  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
1847 
1848 
1849  // Find the clusters that originally had the cells
1850  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1851  {
1852  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1853 
1854  if(idCell>=0)
1855  {
1856  Int_t idCluster = fOrgClusterCellId[idCell];
1857 
1858  Bool_t set = kTRUE;
1859  for(Int_t icl =0; icl < nClu; icl++)
1860  {
1861  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
1862  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
1863  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
1864  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
1865  }
1866  if( set && idCluster >= 0)
1867  {
1868  clArray.SetAt(idCluster,nClu++);
1869  //printf("******** idCluster %d \n",idCluster);
1870  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
1871 
1872  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
1873 
1874  //Search highest E cluster
1875  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1876  //printf("\t E %f\n",clOrg->E());
1877  if(emax < clOrg->E())
1878  {
1879  emax = clOrg->E();
1880  idMax = idCluster;
1881  }
1882  }
1883  }
1884  }// cell loop
1885 
1886  if(nClu==0 || nLabTotOrg == 0)
1887  {
1888  //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());
1889  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
1890  //printf("\n");
1891  }
1892 
1893  // Put the first in the list the cluster with highest energy
1894  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1895  {
1896  Int_t maxIndex = -1;
1897  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1898  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1899  {
1900  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1901  }
1902 
1903  if(firstCluster >=0 && idMax >=0)
1904  {
1905  clArray.SetAt(idMax,0);
1906  clArray.SetAt(firstCluster,maxIndex);
1907  }
1908  }
1909 
1910  // Get the labels list in the original clusters, assign all to the new cluster
1911  TArrayI clMCArray(nLabTotOrg) ;
1912  clMCArray.Reset();
1913 
1914  Int_t nLabTot = 0;
1915  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1916  {
1917  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
1918  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
1919  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1920  Int_t nLab = clOrg->GetNLabels();
1921 
1922  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
1923  {
1924  Int_t lab = clOrg->GetLabelAt(iLab) ;
1925  if(lab>=0)
1926  {
1927  Bool_t set = kTRUE;
1928  //printf("\t \t Set Label %d \n", lab);
1929  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1930  {
1931  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
1932  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
1933  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
1934  }
1935  if( set ) clMCArray.SetAt(lab,nLabTot++);
1936  }
1937  }
1938  }// cluster loop
1939 
1940  // Set the final list of labels
1941 
1942  clus->SetLabel(clMCArray.GetArray(), nLabTot);
1943 
1944 // printf("Final list of labels for new cluster : \n");
1945 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
1946 // {
1947 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
1948 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
1949 // printf(" org %d ",label);
1950 // RemapMCLabelForAODs(label);
1951 // printf(" new %d \n",label);
1952 // }
1953 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
1954 }
1955 
1956 //____________________________________________________________
1957 // Init geometry, create list of output clusters.
1958 //____________________________________________________________
1960 {
1961  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1962 
1963  if(fOutputAODBranchName.Length()==0)
1964  {
1965  fOutputAODBranchName = "newEMCALClustersArray";
1966  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
1967  }
1968 
1970 
1971  if( fFillAODFile )
1972  {
1973  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1974 
1975  //fOutputAODBranch->SetOwner(kFALSE);
1976 
1977  AddAODBranch("TClonesArray", &fOutputAODBranch);
1978  }
1979 }
1980 
1981 //_______________________________________________________
1992 //_______________________________________________________
1994 {
1995  //-------
1996  // Step 1
1997 
1998  // Remove the contents of AOD branch output list set in the previous event
1999  fOutputAODBranch->Clear("C");
2000 
2001  LoadBranches();
2002 
2003  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
2004  // If we need a centrality bin, we select only those events in the corresponding bin.
2005  if( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
2006  {
2007  Float_t cen = GetEventCentrality();
2008  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
2009  }
2010 
2011  // Intermediate array with new clusters : init the array only once or clear from previous event
2012  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
2013  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
2014 
2015  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
2016 
2017  // Get the event, do some checks and settings
2018  CheckAndGetEvent() ;
2019 
2020  if (!fEvent)
2021  {
2022  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
2023  return ;
2024  }
2025 
2026  // Init pointers, geometry, clusterizer, ocdb, aodb
2027 
2028  if(fAccessOCDB) AccessOCDB();
2029  if(fAccessOADB) AccessOADB(); // only once
2030 
2032 
2033  //-------
2034  // Step 2
2035 
2036  // Make clusters
2038  else ClusterizeCells() ;
2039 
2040  //-------
2041  // Step 3
2042 
2044 }
2045 
2046 
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
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.