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