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