AliPhysics  6bc8652 (6bc8652)
 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  // Only 1 L1 phase correction possible, except special cases
434  TString pass2 = "pass1";
435 
436  if ( pass=="muon_calo_pass1" && runnumber > 209121 && runnumber < 244284 )
437  pass2 = "pass0"; // period LHC15a-m
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  nClustersOrg++;
646 
647  if(!clus->IsEMCAL()) continue;
648 
649  Int_t label = clus->GetLabel();
650  Int_t label2 = -1 ;
651  if (clus->GetNLabels() >=2 ) label2 = clus->GetLabelAt(1) ;
652 
653  //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",
654  // i, clus->GetID(), clus->E(), clus->GetTOF()*1e9, clus->GetNCells(), clus->GetNLabels(), label);
655  //
656  //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++)
657  // printf("%d) Label %d, E dep frac %0.2f; ",
658  // imc, clus->GetLabelAt(imc),clus->GetClusterMCEdepFraction(imc));
659  //if(clus->GetNLabels() > 0) printf("\n");
660 
661  UShort_t * index = clus->GetCellsAbsId() ;
662  for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
663  {
664  fOrgClusterCellId[index[icell]] = i;
665  fCellTime[index[icell]] = clus->GetTOF();
666  fCellMatchdEta[index[icell]] = clus->GetTrackDz();
667  fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
668 
670  {
671  fCellLabels[index[icell]] = label;
672  fCellSecondLabels[index[icell]] = label2;
673  }
674  }
675  }
676  }
677 
678  // Transform CaloCells into Digits
679 
680  Int_t idigit = 0;
681  Int_t id = -1;
682  Float_t amp = -1;
683  Double_t time = -1;
684 
685  AliVCaloCells * cells = fEvent->GetEMCALCells();
686 
687  Int_t bc = InputEvent()->GetBunchCrossNumber();
688 
689  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
690  {
691  // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
692  id = cells->GetCellNumber(icell);
693  Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
694  time-=fConstantTimeShift*1e-9; // only in case of simulations done before 2015
695 
696  // Do not include cells with too low energy, nor exotic cell
697  if( amp < fRecParam->GetMinECut() ||
698  time > fRecParam->GetTimeMax() ||
699  time < fRecParam->GetTimeMin() ) accept = kFALSE;
700 
701  // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
703  {
704  time = fCellTime[id];
705  //printf("cell %d time org %f - ",id, time*1.e9);
706  fRecoUtils->RecalibrateCellTime(id,bc,time);
707  //printf("recal %f\n",time*1.e9);
708  }
709 
710  //Exotic?
711  if (accept && fRecoUtils->IsExoticCell(id,cells,bc))
712  accept = kFALSE;
713 
714  if( !accept )
715  {
716  AliDebug(2,Form("Remove channel absId %d, index %d of %d, amp %f, time %f",
717  id,icell, cells->GetNumberOfCells(), amp, time*1.e9));
718  continue;
719  }
720 
721  //
722  // MC
723  //
724 
725  // Old way
726  Int_t mcLabel = cells->GetMCLabel(icell);
727  Float_t eDep = amp;
728 
729  //printf("--- Selected cell %d--- amp %f\n",id,amp);
730 
731  // New way
732  TArrayI labeArr(0);
733  TArrayF eDepArr(0);
734  Int_t nLabels = 0;
735 
737  {
738  // Old way to recover/set the cell MC label
739  // Only possibility for old Run1 productions
740 
741  if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
742 
743  else if( fSetCellMCLabelFromCluster == 0 &&
745 
746  else mcLabel = -1; // found later
747 
748  // Last parameter of the digit object should be MC deposited energy,
749  // since it is not available in aliroot before year 2016, add just the cell amplitude so that
750  // we give more weight to the MC label of the cell with highest energy in the cluster
751 
752  Float_t efrac = cells->GetEFraction(icell);
753 
754  // When checking the MC of digits, give weight to cells with embedded signal
755  if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
756 
757  eDep *= efrac ;
758  }
759  else // fSetCellMCLabelFromEdepFrac = true
760  {
761  // New way, valid only for MC productions with aliroot > v5-07-21
762  mcLabel = -1;
763 
764  // Map the digit to cell index for later to calculate the cell MC energy deposition map
765  fCellLabels[id] = idigit;
766  //printf("\t absId %d, idigit %d\n",id,idigit);
767 
768  if(fOrgClusterCellId[id] < 0) continue; // index can be negative if noisy cell that did not form cluster
769 
770  AliVCluster *clus = 0;
771  Int_t iclus = fOrgClusterCellId[id];
772 
773  if(iclus < 0)
774  {
775  AliInfo("Negative original cluster index, skip \n");
776  continue;
777  }
778 
779  if(aodIH && aodIH->GetEventToMerge()) //Embedding
780  clus = aodIH->GetEventToMerge()->GetCaloCluster(iclus); //Get clusters directly from embedded signal
781  else
782  clus = fEvent->GetCaloCluster(iclus);
783 
784  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(id, clus, MCEvent(), amp, labeArr, eDepArr);
785  nLabels = labeArr.GetSize();
786 
787  } // cell MC label, new
788 
789  //
790  // Create the digit
791  //
792  if(amp <= 0.01) continue ; // accept if > 10 MeV
793 
794  //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);
795 
796  AliEMCALDigit* digit = new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, eDep);
797 
798  if(nLabels > 0)
799  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
800 
801  idigit++;
802 
803  }
804 
805  fDigitsArr->Sort();
806 
807  //-------------------------------------------------------------------------------------
808  // Do the clusterization
809  //-------------------------------------------------------------------------------------
810 
811  fClusterizer->Digits2Clusters("");
812 
813  //-------------------------------------------------------------------------------------
814  // Transform the recpoints into AliVClusters
815  //-------------------------------------------------------------------------------------
816 
818 
819  if(!fCaloClusterArr)
820  {
821  AliWarning(Form("No array with CaloClusters, input RecPoints entries %d",fClusterArr->GetEntriesFast()));
822  return;
823  }
824 
825  AliDebug(1,Form("N clusters: before recluster %d, after recluster %d, recpoints %d",
826  nClustersOrg, fCaloClusterArr->GetEntriesFast(),fClusterArr->GetEntriesFast()));
827 
828 // if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
829 // {
830 // AliInfo("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
831 // fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
832 // fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
833 // }
834 }
835 
836 
837 //_____________________________________________________
839 //_____________________________________________________
841 {
842  AliVCaloCells *cells = fEvent->GetEMCALCells();
843  Double_t cellAmplitude = 0;
844  Double_t cellTime = 0;
845  Short_t cellNumber = 0;
846  Int_t cellMCLabel = 0;
847  Double_t cellEFrac = 0;
848  Int_t nClustersOrg = 0;
849 
850  // Fill the array with the EMCAL clusters, copy them
851  for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
852  {
853  AliVCluster *clus = fEvent->GetCaloCluster(i);
854  if(clus->IsEMCAL())
855  {
856  //recalibrate/remove bad channels/etc if requested
857  if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
858  {
859  continue;
860  }
861 
862  if(fRecoUtils->IsRecalibrationOn())
863  {
864  //Calibrate cluster
865  fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
866 
867  //CalibrateCells
868  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
869  {
870  if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
871  break;
872 
873  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
874  fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
875  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
876 
877  //Do not include bad channels found in analysis?
878  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
879  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
880  continue;
881 
882  cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
883 
884  }// cells loop
885  }// recalibrate
886 
887  //Cast to ESD or AOD, needed to create the cluster array
888  AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
889  AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
890 
891  if (esdCluster)
892  {
893  fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
894  }//ESD
895  else if(aodCluster)
896  {
897  fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
898  }//AOD
899  else
900  AliWarning("Wrong CaloCluster type?");
901 
902  nClustersOrg++;
903  }
904  }
905 
906  //Do the unfolding
907  fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
908 
909  //CLEAN-UP
910  fUnfolder->Clear();
911 }
912 
913 //_____________________________________________________
915 //_____________________________________________________
917 {
918  AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
919  Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
920 
921  AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
922  aodEMcells.CreateContainer(nEMcell);
923  aodEMcells.SetType(AliVCaloCells::kEMCALCell);
924  Double_t calibFactor = 1.;
925  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
926  {
927  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
928  fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
929  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
930 
931  if(fRecoUtils->IsRecalibrationOn())
932  {
933  calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
934  }
935 
936  if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
937  { //Channel is not declared as bad
938  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
939  eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
940  }
941  else
942  {
943  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
944  }
945  }
946 
947  aodEMcells.Sort();
948 }
949 
950 //__________________________________________________
952 //__________________________________________________
954 {
955  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
956  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
957 
958  Double_t pos[3] ;
959  Double_t covVtx[6];
960  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
961 
962  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
963  if(!header) AliFatal("Not a standard AOD");
964  header->SetRunNumber(fEvent->GetRunNumber());
965 
966  if(esdevent)
967  {
968  TTree* tree = fInputHandler->GetTree();
969  if (tree)
970  {
971  TFile* file = tree->GetCurrentFile();
972  if (file) header->SetESDFileName(file->GetName());
973  }
974  }
975  else if (aodevent) {
976  AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
977  if(!aodheader) AliFatal("Not a standard AOD");
978  header->SetESDFileName(aodheader->GetESDFileName());
979  }
980 
981  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
982  header->SetOrbitNumber(fEvent->GetOrbitNumber());
983  header->SetPeriodNumber(fEvent->GetPeriodNumber());
984  header->SetEventType(fEvent->GetEventType());
985 
986  // Centrality
988  {
989  if(fEvent->GetCentrality())
990  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
991  else
992  header->SetCentrality(0);
993  }
994 
995  // Trigger
996  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
997 
998  header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
999 
1000  header->SetTriggerMask(fEvent->GetTriggerMask());
1001  header->SetTriggerCluster(fEvent->GetTriggerCluster());
1002 
1003  if (esdevent)
1004  {
1005  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
1006  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
1007  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
1008  }
1009  else if (aodevent)
1010  {
1011  header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
1012  header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
1013  header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
1014  }
1015 
1016  header->SetMagneticField(fEvent->GetMagneticField());
1017  //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
1018 
1019  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
1020  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
1021  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
1022  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
1023  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
1024 
1025  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
1026  Float_t diamcov[3];
1027  fEvent->GetDiamondCovXY(diamcov);
1028  header->SetDiamond(diamxy,diamcov);
1029  if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
1030  else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
1031 
1032  //
1033  Int_t nVertices = 1 ;/* = prim. vtx*/;
1034  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1035 
1036  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1037 
1038  // Access to the AOD container of vertices
1039  TClonesArray &vertices = *(AODEvent()->GetVertices());
1040  Int_t jVertices=0;
1041 
1042  // Add primary vertex. The primary tracks will be defined
1043  // after the loops on the composite objects (V0, cascades, kinks)
1044  fEvent->GetPrimaryVertex()->GetXYZ(pos);
1045  Float_t chi = 0;
1046  if (esdevent)
1047  {
1048  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1049  chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
1050  }
1051  else if (aodevent)
1052  {
1053  aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1054  chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
1055  }
1056 
1057  AliAODVertex * primary = new(vertices[jVertices++])
1058  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
1059  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
1060  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
1061 }
1062 
1063 //___________________________________________________________
1067 //___________________________________________________________
1069 {
1070  // First recalculate track-matching for the new clusters
1071  if(fDoTrackMatching)
1072  {
1073  fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom,MCEvent());
1074  }
1075 
1076  // Put the new clusters in the AOD list
1077 
1078  Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
1079 
1080  for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
1081  {
1082  AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1083 
1084  newCluster->SetID(i);
1085 
1086  // Correct cluster energy non linearity
1087  newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
1088 
1089  //Add matched track
1090  if(fDoTrackMatching)
1091  {
1092  Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1093  if(trackIndex >= 0)
1094  {
1095  newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1096  AliDebug(2,Form("Matched Track index %d to new cluster %d",trackIndex,i));
1097  }
1098 
1099  Float_t dR = 999., dZ = 999.;
1100  fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dZ,dR);
1101  newCluster->SetTrackDistance(dR,dZ);
1102  }
1103  else
1104  {// Assign previously assigned matched track in reco, very very rough
1105  Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1106  newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1107  }
1108 
1109  //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1110  //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1111  //printf("\n");
1112  //
1113  //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",
1114  // i, newCluster->GetID(), newCluster->E(), newCluster->GetTOF()*1e9, newCluster->GetNCells(), newCluster->GetNLabels(), newCluster->GetLabel());
1115  //
1116  //for(Int_t imc = 0; imc < newCluster->GetNLabels(); imc++)
1117  // printf("%d) Label %d, E dep frac %0.2f; ",
1118  // imc,newCluster->GetLabelAt(imc),newCluster->GetClusterMCEdepFraction(imc));
1119  //if(newCluster->GetNLabels() > 0) printf("\n");
1120 
1121  // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1122  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
1123 
1124  new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
1125 
1126  AliDebug(2,Form("New cluster %d of %d, energy %f, mc label %d",
1127  newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()));
1128 
1129  } // cluster loop
1130 
1131  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1132 }
1133 
1134 
1135 //________________________________________________________________
1137 //________________________________________________________________
1139 {
1140  if(fUseAliCentrality)
1141  {
1142  if(GetCentrality()) return GetCentrality()->GetCentralityPercentile(fCentralityClass) ;
1143  else return -1. ;
1144  }
1145  else
1146  {
1147  if(GetMultSelCen()) return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass,kTRUE) ;
1148  else return -1. ;
1149  }
1150 
1151  return -1.;
1152 }
1153 
1154 //_______________________________________________
1156 //_______________________________________________
1158 {
1159  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1160  {
1161  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null");
1162  return TString("");
1163  }
1164 
1165  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1166  {
1167  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null");
1168  return TString("");
1169  }
1170 
1171  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1172  if (pass.Contains("ass1")) return TString("pass1");
1173  else if (pass.Contains("ass2")) return TString("pass2");
1174  else if (pass.Contains("ass3")) return TString("pass3");
1175  else if (pass.Contains("ass4")) return TString("pass4");
1176  else if (pass.Contains("ass5")) return TString("pass5");
1177  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1178  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1179  {
1180  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1181  return TString("pass1");
1182  }
1183  else if (pass.Contains("LHC14a1a"))
1184  {
1185  AliInfo("Enable EMCal energy calibration for this MC production!!");
1186 
1187  fRecoUtils->SwitchOnRecalibration();
1188 
1189  return TString("LHC14a1a");
1190  }
1191 
1192  // No condition fullfilled, give a default value
1193  AliInfo("Pass number string not found");
1194 
1195  return TString("");
1196 }
1197 
1198 //_________________________________________
1201 //_________________________________________
1203 {
1204  if(fDebug >=0) (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1205 
1206  fOADBSet = kFALSE;
1207  if(fOADBFilePath == "") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
1208 
1209  fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1210 
1211  if(!fRecParam) fRecParam = new AliEMCALRecParam;
1212  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1213 
1214  if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1215  if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1216  if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1217  fRejectBelowThreshold = kFALSE;
1218 
1219  // Centrality
1220  if(fCentralityClass == "") fCentralityClass = "V0M";
1221 
1222  if (fOCDBpath == "") fOCDBpath = "raw://" ;
1223  if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1224 
1225  if(gROOT->LoadMacro(fConfigName) >=0)
1226  {
1227  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
1228  AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1229  fGeomName = clus->fGeomName;
1231  fOCDBpath = clus->fOCDBpath;
1232  fAccessOCDB = clus->fAccessOCDB;
1233  fRecParam = clus->fRecParam;
1234  fJustUnfold = clus->fJustUnfold;
1235  fFillAODFile = clus->fFillAODFile;
1236  fRecoUtils = clus->fRecoUtils;
1237  fConfigName = clus->fConfigName;
1238  fMaxEvent = clus->fMaxEvent;
1241  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1243  fCentralityBin[0] = clus->fCentralityBin[0];
1244  fCentralityBin[1] = clus->fCentralityBin[1];
1246  }
1247 }
1248 
1249 //_______________________________________________________
1252 //_______________________________________________________
1254 {
1255  if (fJustUnfold)
1256  {
1257  // init the unfolding afterburner
1258  delete fUnfolder;
1259  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1260  return;
1261  }
1262 
1263  //First init the clusterizer
1264  delete fClusterizer;
1265  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1266  fClusterizer = new AliEMCALClusterizerv1 (fGeom);
1267  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1268  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1269  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1270  {
1271  fClusterizer = new AliEMCALClusterizerNxN(fGeom);
1272  fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1273  fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1274  }
1275  else
1276  {
1277  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1278  }
1279 
1280  //Now set the parameters
1281  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1282  fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1283  fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1284  fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1285  fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1286  fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1287  fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1288  fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1289  fClusterizer->SetTimeCalibration ( fRecParam->IsTimeCalibrationOn() );
1290  fClusterizer->SetInputCalibrated ( kTRUE );
1291  fClusterizer->SetJustClusters ( kTRUE );
1292 
1293  // Initialize the cluster rec points and digits arrays and get them.
1294  fClusterizer->SetOutput(0);
1295  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1296  fDigitsArr = fClusterizer->GetDigits();
1297 
1298  //In case of unfolding after clusterization is requested, set the corresponding parameters
1299  if(fRecParam->GetUnfold())
1300  {
1301  Int_t i=0;
1302  for (i = 0; i < 8; i++)
1303  {
1304  fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1305  }//end of loop over parameters
1306 
1307  for (i = 0; i < 3; i++)
1308  {
1309  fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1310  fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1311  }//end of loop over parameters
1312  fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1313  fClusterizer->InitClusterUnfolding();
1314 
1315  }// to unfold
1316 }
1317 
1318 //_________________________________________________
1322 //_________________________________________________
1324 {
1325  if(fGeomMatrixSet) return;
1326 
1327  Int_t runnumber = InputEvent()->GetRunNumber() ;
1328  if (!fGeom)
1329  {
1330  if(fGeomName=="")
1331  {
1332  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(runnumber);
1333  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fGeom->GetName(),runnumber));
1334  }
1335  else
1336  {
1337  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1338  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeomName.Data()));
1339  }
1340 
1341  // Init geometry, I do not like much to do it like this ...
1342  if(fImportGeometryFromFile && !gGeoManager)
1343  {
1344  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1345  {
1346  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1347  if (runnumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1348  else if (runnumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1349  else if (runnumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1350  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >=2015
1351  }
1352 
1353  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1354 
1355  TGeoManager::Import(fImportGeometryFilePath) ;
1356  }
1357 
1358  AliDebug(1,Form("Init for run=%d",runnumber));
1359  if (!gGeoManager) AliDebug(1,"Careful!, gGeoManager not loaded, load misalign matrices");
1360  } // geometry pointer did not exist before
1361 
1362  if(fLoadGeomMatrices)
1363  {
1364  AliInfo("Load user defined EMCAL geometry matrices");
1365 
1366  // OADB if available
1367  AliOADBContainer emcGeoMat("AliEMCALgeo");
1368  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1369  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1370 
1371  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1372  {
1373 
1374  if (!fGeomMatrix[mod]) // Get it from OADB
1375  {
1376  AliDebug(2,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
1377  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1378 
1379  fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1380  }
1381 
1382  if(fGeomMatrix[mod])
1383  {
1384  if(DebugLevel() > 1)
1385  fGeomMatrix[mod]->Print();
1386 
1387  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1388  }
1389  else if(gGeoManager)
1390  {
1391  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
1392  fGeom->SetMisalMatrix(fGeom->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
1393  }
1394  else
1395  {
1396  AliError(Form("Alignment atrix for SM %d is not available",mod));
1397  }
1398 
1399  fGeomMatrixSet=kTRUE;
1400 
1401  }//SM loop
1402  }//Load matrices
1403  else if(!gGeoManager)
1404  {
1405  AliInfo("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1406  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1407  if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1408  {
1409  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
1410  }//AOD
1411  else
1412  {
1413  AliDebug(1,"Load Misaligned matrices");
1414 
1415  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1416 
1417  if(!esd)
1418  {
1419  AliError("This event does not contain ESDs?");
1420  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1421  return;
1422  }
1423 
1424  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1425  {
1426  if(DebugLevel() > 1)
1427  esd->GetEMCALMatrix(mod)->Print();
1428 
1429  if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1430 
1431  }
1432 
1433  fGeomMatrixSet=kTRUE;
1434 
1435  } // ESD
1436  } // Load matrices from Data
1437 }
1438 
1439 //____________________________________________________
1444 //____________________________________________________
1446 {
1447  if(!fRemoveExoticEvents) return kFALSE;
1448 
1449  // Loop on cells
1450 
1451  AliVCaloCells * cells = fEvent->GetEMCALCells();
1452  Float_t totCellE = 0;
1453  Int_t bc = InputEvent()->GetBunchCrossNumber();
1454 
1455  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1456  {
1457  Float_t ecell = 0 ;
1458  Double_t tcell = 0 ;
1459 
1460  Int_t absID = cells->GetCellNumber(icell);
1461  Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1462  tcell-=fConstantTimeShift*1e-9;// Only for MC simulations done before 2015
1463 
1464  if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1465  }
1466 
1467  // TString triggerclasses = event->GetFiredTriggerClasses();
1468  // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1469  // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1470  // return;
1471  //
1472 
1473  //printf("TotE cell %f\n",totCellE);
1474  if(totCellE < 1) return kTRUE;
1475 
1476  return kFALSE;
1477 }
1478 
1479 //________________________________________________________________
1481 //________________________________________________________________
1483 {
1484  if(!fRemoveLEDEvents) return kFALSE;
1485 
1486  //check events of LHC11a period
1487  if(run < 146858 || run > 146860) return kFALSE ;
1488 
1489  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1490  Int_t ncellsSM3 = 0;
1491  AliVCaloCells * cells = fEvent->GetEMCALCells();
1492  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1493  {
1494  if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1495  }
1496 
1497  TString triggerclasses = fEvent->GetFiredTriggerClasses();
1498 
1499  Int_t ncellcut = 21;
1500  if(triggerclasses.Contains("EMC")) ncellcut = 35;
1501 
1502  if( ncellsSM3 >= ncellcut)
1503  {
1504  AliInfo(Form("Reject event %d with ncells in SM3 %d",(Int_t)Entry(),ncellsSM3));
1505  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1506  return kTRUE;
1507  }
1508 
1509  return kFALSE;
1510 }
1511 
1512 //_______________________________________________________
1516 //_______________________________________________________
1518 {
1519  Int_t j = 0;
1520  for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1521  {
1522  AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1523 
1524  const Int_t ncells = recPoint->GetMultiplicity();
1525  Int_t ncellsTrue = 0;
1526 
1527  if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1528 
1529  // cells and their amplitude fractions
1530  UShort_t absIds[ncells];
1531  Double32_t ratios[ncells];
1532 
1533  //For later check embedding
1534  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1535 
1536  Float_t clusterE = 0;
1537  for (Int_t c = 0; c < ncells; c++)
1538  {
1539  AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1540 
1541  absIds[ncellsTrue] = digit->GetId();
1542  ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1543 
1544  // In case of unfolding, remove digits with unfolded energy too low
1545  if(fSelectCell)
1546  {
1547  if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1548  {
1549 
1550  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
1551  recPoint->GetEnergiesList()[c],digit->GetAmplitude()));
1552  continue;
1553  } // if cuts
1554  }// Select cells
1555 
1556  //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1557  clusterE +=recPoint->GetEnergiesList()[c];
1558 
1559  // In case of embedding, fill ratio with amount of signal,
1560  if (aodIH && aodIH->GetMergeEvents())
1561  {
1562  //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1563  AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1564  AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1565 
1566  Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1567  //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1568  Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1569  //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1570 
1571  if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1572  //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1573 
1574  }//Embedding
1575 
1576  ncellsTrue++;
1577 
1578  }// cluster cell loop
1579 
1580  if (ncellsTrue < 1)
1581  {
1582  AliDebug(2,Form("Skipping cluster with no cells avobe threshold E = %f, ncells %d",
1583  recPoint->GetEnergy(), ncells));
1584  continue;
1585  }
1586 
1587  //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1588 
1589  if(clusterE < fRecParam->GetClusteringThreshold())
1590  {
1591  AliDebug(2,Form("Remove cluster with energy below seed threshold %f",clusterE));
1592  continue;
1593  }
1594 
1595  TVector3 gpos;
1596  Float_t g[3];
1597 
1598  // calculate new cluster position
1599 
1600  recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1601  recPoint->GetGlobalPosition(gpos);
1602  gpos.GetXYZ(g);
1603 
1604  // create a new cluster
1605 
1606  (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1607  AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1608  j++;
1609  clus->SetType(AliVCluster::kEMCALClusterv1);
1610  clus->SetE(clusterE);
1611  clus->SetPosition(g);
1612  clus->SetNCells(ncellsTrue);
1613  clus->SetCellsAbsId(absIds);
1614  clus->SetCellsAmplitudeFraction(ratios);
1615  clus->SetChi2(-1); //not yet implemented
1616  clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1617  clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1618  clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1619 
1620  if(ncells == ncellsTrue)
1621  {
1622  Float_t elipAxis[2];
1623  recPoint->GetElipsAxis(elipAxis);
1624  clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1625  clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1626  clus->SetDispersion(recPoint->GetDispersion());
1627  }
1628  else if(fSelectCell)
1629  {
1630  // In case some cells rejected, in unfolding case, recalculate
1631  // shower shape parameters and position
1632  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
1633 
1634  AliVCaloCells* cells = 0x0;
1635  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1636  else cells = InputEvent()->GetEMCALCells();
1637 
1638  fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1639  fRecoUtils->RecalculateClusterPID(clus);
1640  fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1641  }
1642 
1643  // MC
1644 
1647  else
1648  {
1649  //
1650  // Normal case, trust what the clusterizer has found
1651  //
1652  Int_t parentMult = 0;
1653  Int_t *parentList = recPoint->GetParents(parentMult);
1654  Float_t *parentListDE = recPoint->GetParentsDE(); // deposited energy
1655 
1656  clus->SetLabel(parentList, parentMult);
1658  clus->SetClusterMCEdepFractionFromEdepArray(parentListDE);
1659 
1660  //
1661  // Set the cell energy deposition fraction map:
1662  //
1663  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
1664  {
1665  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
1666 
1667  // Get the digit that originated this cell cluster
1668  AliVCaloCells* cells = 0x0;
1669  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1670  else cells = InputEvent()->GetEMCALCells();
1671 
1672  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
1673  {
1674  Int_t idigit = fCellLabels[absIds[icell]];
1675 
1676  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
1677 
1678  // Find the 4 MC labels that contributed to the cluster and their
1679  // deposited energy in the current digit
1680 
1681  mcEdepFracPerCell[icell] = 0; // init
1682 
1683  Int_t nparents = dig->GetNiparent();
1684  if ( nparents > 0 )
1685  {
1686  Int_t digLabel =-1 ;
1687  Float_t edep = 0 ;
1688  Float_t edepTot = 0 ;
1689  Float_t mcEDepFrac[4] = {0,0,0,0};
1690 
1691  // all parents in digit
1692  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
1693  {
1694  digLabel = dig->GetIparent (jndex+1);
1695  edep = dig->GetDEParent(jndex+1);
1696  edepTot += edep;
1697 
1698  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
1699  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
1700  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
1701  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
1702  } // all prarents in digit
1703 
1704  // Divide energy deposit by total deposited energy
1705  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
1706  if(edepTot > 0.01)
1707  {
1708  mcEdepFracPerCell[icell] = clus->PackMCEdepFraction(mcEDepFrac);
1709  }
1710  } // at least one parent label in digit
1711  } // cell in cluster loop
1712 
1713  clus->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
1714 
1715  delete [] mcEdepFracPerCell;
1716 
1717  } // at least one parent in cluster, do the cell primary packing
1718  }
1719  } // recPoints loop
1720 }
1721 
1722 //_____________________________________________________________________
1725 //_____________________________________________________________________
1727 {
1728  if(label < 0) return ;
1729 
1730  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
1731  if(!evt) return ;
1732 
1733  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1734  if(!arr) return ;
1735 
1736  if(label < arr->GetEntriesFast())
1737  {
1738  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1739  if(!particle) return ;
1740 
1741  if(label == particle->Label()) return ; // label already OK
1742  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
1743  }
1744  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
1745 
1746  // loop on the particles list and check if there is one with the same label
1747  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
1748  {
1749  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1750  if(!particle) continue ;
1751 
1752  if(label == particle->Label())
1753  {
1754  label = ind;
1755  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
1756  return;
1757  }
1758  }
1759 
1760  label = -1;
1761 
1762  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
1763 }
1764 
1765 //________________________________________________
1767 //________________________________________________
1769 {
1770  for(Int_t j = 0; j < fgkNEMCalCells; j++)
1771  {
1772  fOrgClusterCellId[j] = -1;
1773  fCellLabels[j] = -1;
1774  fCellSecondLabels[j] = -1;
1775  fCellTime[j] = 0.;
1776  fCellMatchdEta[j] = -999;
1777  fCellMatchdPhi[j] = -999;
1778  }
1779 }
1780 
1781 //_____________________________________________________________________________________________________
1785 //_____________________________________________________________________________________________________
1787  AliAODCaloCluster * clus)
1788 {
1789  Int_t parentMult = 0;
1790  Int_t *parentList = recPoint->GetParents(parentMult);
1791  clus->SetLabel(parentList, parentMult);
1792 
1793  //Write the second major contributor to each MC cluster.
1794  Int_t iNewLabel ;
1795  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1796  {
1797 
1798  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1799  if(idCell>=0)
1800  {
1801  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
1802  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1803  {
1804  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
1805  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
1806  }
1807  if (iNewLabel == 1)
1808  {
1809  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1810  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1811  {
1812  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1813  }
1814 
1815  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1816  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1817  delete [] newLabelArray;
1818  }
1819  } // positive cell number
1820  }
1821 }
1822 
1823 //___________________________________________________________________________________________________
1827 //___________________________________________________________________________________________________
1829 {
1830  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1831  clArray.Reset();
1832  Int_t nClu = 0;
1833  Int_t nLabTotOrg = 0;
1834  Float_t emax = -1;
1835  Int_t idMax = -1;
1836 
1837  AliVEvent * event = fEvent;
1838 
1839  // In case of embedding the MC signal is in the added event
1840  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1841  if(aodIH && aodIH->GetEventToMerge()) //Embedding
1842  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
1843 
1844 
1845  // Find the clusters that originally had the cells
1846  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1847  {
1848  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1849 
1850  if(idCell>=0)
1851  {
1852  Int_t idCluster = fOrgClusterCellId[idCell];
1853 
1854  Bool_t set = kTRUE;
1855  for(Int_t icl =0; icl < nClu; icl++)
1856  {
1857  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
1858  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
1859  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
1860  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
1861  }
1862  if( set && idCluster >= 0)
1863  {
1864  clArray.SetAt(idCluster,nClu++);
1865  //printf("******** idCluster %d \n",idCluster);
1866  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
1867 
1868  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
1869 
1870  //Search highest E cluster
1871  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1872  //printf("\t E %f\n",clOrg->E());
1873  if(emax < clOrg->E())
1874  {
1875  emax = clOrg->E();
1876  idMax = idCluster;
1877  }
1878  }
1879  }
1880  }// cell loop
1881 
1882  if(nClu==0 || nLabTotOrg == 0)
1883  {
1884  //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());
1885  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
1886  //printf("\n");
1887  }
1888 
1889  // Put the first in the list the cluster with highest energy
1890  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1891  {
1892  Int_t maxIndex = -1;
1893  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1894  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1895  {
1896  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1897  }
1898 
1899  if(firstCluster >=0 && idMax >=0)
1900  {
1901  clArray.SetAt(idMax,0);
1902  clArray.SetAt(firstCluster,maxIndex);
1903  }
1904  }
1905 
1906  // Get the labels list in the original clusters, assign all to the new cluster
1907  TArrayI clMCArray(nLabTotOrg) ;
1908  clMCArray.Reset();
1909 
1910  Int_t nLabTot = 0;
1911  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1912  {
1913  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
1914  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
1915  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1916  Int_t nLab = clOrg->GetNLabels();
1917 
1918  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
1919  {
1920  Int_t lab = clOrg->GetLabelAt(iLab) ;
1921  if(lab>=0)
1922  {
1923  Bool_t set = kTRUE;
1924  //printf("\t \t Set Label %d \n", lab);
1925  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1926  {
1927  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
1928  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
1929  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
1930  }
1931  if( set ) clMCArray.SetAt(lab,nLabTot++);
1932  }
1933  }
1934  }// cluster loop
1935 
1936  // Set the final list of labels
1937 
1938  clus->SetLabel(clMCArray.GetArray(), nLabTot);
1939 
1940 // printf("Final list of labels for new cluster : \n");
1941 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
1942 // {
1943 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
1944 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
1945 // printf(" org %d ",label);
1946 // RemapMCLabelForAODs(label);
1947 // printf(" new %d \n",label);
1948 // }
1949 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
1950 }
1951 
1952 //____________________________________________________________
1953 // Init geometry, create list of output clusters.
1954 //____________________________________________________________
1956 {
1957  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1958 
1959  if(fOutputAODBranchName.Length()==0)
1960  {
1961  fOutputAODBranchName = "newEMCALClustersArray";
1962  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
1963  }
1964 
1966 
1967  if( fFillAODFile )
1968  {
1969  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1970 
1971  //fOutputAODBranch->SetOwner(kFALSE);
1972 
1973  AddAODBranch("TClonesArray", &fOutputAODBranch);
1974  }
1975 }
1976 
1977 //_______________________________________________________
1988 //_______________________________________________________
1990 {
1991  //-------
1992  // Step 1
1993 
1994  // Remove the contents of AOD branch output list set in the previous event
1995  fOutputAODBranch->Clear("C");
1996 
1997  LoadBranches();
1998 
1999  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
2000  // If we need a centrality bin, we select only those events in the corresponding bin.
2001  if( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
2002  {
2003  Float_t cen = GetEventCentrality();
2004  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
2005  }
2006 
2007  // Intermediate array with new clusters : init the array only once or clear from previous event
2008  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
2009  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
2010 
2011  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
2012 
2013  // Get the event, do some checks and settings
2014  CheckAndGetEvent() ;
2015 
2016  if (!fEvent)
2017  {
2018  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
2019  return ;
2020  }
2021 
2022  // Init pointers, geometry, clusterizer, ocdb, aodb
2023 
2024  if(fAccessOCDB) AccessOCDB();
2025  if(fAccessOADB) AccessOADB(); // only once
2026 
2028 
2029  //-------
2030  // Step 2
2031 
2032  // Make clusters
2034  else ClusterizeCells() ;
2035 
2036  //-------
2037  // Step 3
2038 
2040 }
2041 
2042 
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.