AliPhysics  ebc57ae (ebc57ae)
 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)
69 , fCaloClusterArr(0), fCaloCells(0)
70 , fRecParam(0), fClusterizer(0)
71 , fUnfolder(0), fJustUnfold(kFALSE)
72 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
73 , fFillAODFile(kFALSE), fFillAODHeader(0)
74 , fFillAODCaloCells(0), fRun(-1)
75 , fRecoUtils(0), fConfigName("")
76 , fOrgClusterCellId()
77 , fCellLabels(), fCellSecondLabels(), fCellTime()
78 , fCellMatchdEta(), fCellMatchdPhi()
79 , fRecalibrateWithClusterTime(0)
80 , fMaxEvent(0), fDoTrackMatching(kFALSE), fUpdateCell(0)
81 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
82 , fRejectBelowThreshold(kFALSE)
83 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
84 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
85 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
86 , fConstantTimeShift(0)
87 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
88 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
89 , fSetCellMCLabelFromCluster(0)
90 , fSetCellMCLabelFromEdepFrac(0)
91 , fRemapMCLabelForAODs(0)
92 , fInputFromFilter(0)
93 , fTCardCorrEmulation(0), fTCardCorrClusEnerConserv(0)
94 , fRandom(0)
95 {
96  for(Int_t i = 0; i < 22; i++)
97  {
98  fGeomMatrix[i] = 0;
100  }
101 
102  ResetArrays();
103 
104  fCentralityBin[0] = fCentralityBin[1]=-1;
105 
106  for(Int_t i = 0; i < 4 ; i++)
107  {
108  fTCardCorrInduceEnerFrac [i] = 0 ;
109  fTCardCorrInduceEnerFracP1 [i] = 0 ;
111  }
112 }
113 
114 //______________________________________________________________
116 //______________________________________________________________
118 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
119 , fEvent(0)
120 , fGeom(0), fGeomName("")
121 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
122 , fOCDBpath(""), fAccessOCDB(kFALSE)
123 , fDigitsArr(0), fClusterArr(0)
124 , fCaloClusterArr(0), fCaloCells(0)
125 , fRecParam(0), fClusterizer(0)
126 , fUnfolder(0), fJustUnfold(kFALSE)
127 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
128 , fFillAODFile(kFALSE), fFillAODHeader(0)
129 , fFillAODCaloCells(0), fRun(-1)
130 , fRecoUtils(0), fConfigName("")
131 , fOrgClusterCellId()
132 , fCellLabels(), fCellSecondLabels(), fCellTime()
133 , fCellMatchdEta(), fCellMatchdPhi()
134 , fRecalibrateWithClusterTime(0)
135 , fMaxEvent(0), fDoTrackMatching(kFALSE), fUpdateCell(0)
136 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
137 , fRejectBelowThreshold(kFALSE)
138 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
139 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
140 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
141 , fConstantTimeShift(0)
142 , fCentralityClass(""), fUseAliCentrality(0), fSelectEMCALEvent(0)
143 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
144 , fSetCellMCLabelFromCluster(0)
145 , fSetCellMCLabelFromEdepFrac(0)
146 , fRemapMCLabelForAODs(0)
147 , fInputFromFilter(0)
148 , fTCardCorrEmulation(0), fTCardCorrClusEnerConserv(0)
149 , fRandom(0)
150 {
151  for(Int_t i = 0; i < 22; i++)
152  {
153  fGeomMatrix[i] = 0;
155  }
156 
157  ResetArrays();
158 
159  fCentralityBin[0] = fCentralityBin[1]=-1;
160 
161  for(Int_t i = 0; i < 4 ; i++)
162  {
163  fTCardCorrInduceEnerFrac [i] = 0 ;
164  fTCardCorrInduceEnerFracP1 [i] = 0 ;
166  }
167 }
168 
169 //_______________________________________________________________
171 //_______________________________________________________________
173 {
174  if (fDigitsArr)
175  {
176  fDigitsArr->Clear("C");
177  delete fDigitsArr;
178  }
179 
180  if (fClusterArr)
181  {
182  fClusterArr->Delete();
183  delete fClusterArr;
184  }
185 
186  if (fCaloClusterArr)
187  {
188  fCaloClusterArr->Delete();
189  delete fCaloClusterArr;
190  }
191 
192  if(fClusterizer) delete fClusterizer;
193  if(fUnfolder) delete fUnfolder;
194  if(fRecoUtils) delete fRecoUtils;
195 }
196 
197 //_______________________________________________________
205 //_______________________________________________________________________________
207 {
208  if ( absID < 0 || absID >= 24*48*fGeom->GetNumberOfSuperModules() )
209  return kFALSE;
210 
211  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
212  if (!fGeom->GetCellIndex(absID,imod,iTower,iIphi,iIeta))
213  return kFALSE;
214 
215  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
216 
217  // Do not include bad channels found in analysis,
218  if ( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
219  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi) )
220  return kFALSE;
221 
222  return kTRUE;
223 }
224 
225 //_______________________________________________________
229 //_______________________________________________________
231 {
232  if(!fSelectEMCALEvent) return kTRUE; // accept
233 
234  if(fEMCALEnergyCut <= 0) return kTRUE; // accept
235 
236  Int_t nCluster = fEvent -> GetNumberOfCaloClusters();
237  Int_t bc = fEvent -> GetBunchCrossNumber();
238 
239  for(Int_t icalo = 0; icalo < nCluster; icalo++)
240  {
241  AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
242 
243  if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
244  fRecoUtils->IsGoodCluster(clus,fGeom,fCaloCells,bc))
245  {
246 
247  AliDebug(1, Form("Accept : E %2.2f > %2.2f, nCells %d > %d",
248  clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut));
249 
250  return kTRUE;
251  }
252 
253  }// loop
254 
255  AliDebug(1,"Reject");
256 
257  return kFALSE;
258 }
259 
260 //_______________________________________________
263 //_______________________________________________
265 {
266  //Set it only once
267  if(fOADBSet) return ;
268 
269  Int_t runnumber = InputEvent()->GetRunNumber() ;
270  TString pass = GetPass();
271 
272  AliInfo(Form("Get AODB parameters from EMCAL in %s for run %d, and <%s>",fOADBFilePath.Data(),runnumber,pass.Data()));
273 
274  Int_t nSM = fGeom->GetNumberOfSuperModules();
275 
276  // Bad map
277  if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
278  {
279  AliOADBContainer *contBC=new AliOADBContainer("");
280  contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
281 
282  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
283 
284  if(arrayBC)
285  {
286  fRecoUtils->SwitchOnDistToBadChannelRecalculation();
287  AliInfo("Remove EMCAL bad cells");
288 
289  for (Int_t i=0; i < nSM; ++i)
290  {
291  TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
292 
293  if (hbm)
294  delete hbm;
295 
296  hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
297 
298  if (!hbm)
299  {
300  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
301  continue;
302  }
303 
304  hbm->SetDirectory(0);
305  fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
306 
307  } // loop
308  } else AliInfo("Do NOT remove EMCAL bad channels"); // run array
309 
310  delete contBC;
311  } // Remove bad
312 
313  // Energy Recalibration
314  if(fRecoUtils->IsRecalibrationOn())
315  {
316  AliOADBContainer *contRF=new AliOADBContainer("");
317 
318  contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
319 
320  TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
321 
322  if(recal)
323  {
324  TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
325 
326  if(recalpass)
327  {
328  TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
329 
330  if(recalib)
331  {
332  AliInfo("Recalibrate EMCAL");
333  for (Int_t i=0; i<nSM; ++i)
334  {
335  TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
336 
337  if (h)
338  delete h;
339 
340  h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
341 
342  if (!h)
343  {
344  AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
345  continue;
346  }
347 
348  h->SetDirectory(0);
349 
350  fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
351  } // SM loop
352  } else AliInfo("Do NOT recalibrate EMCAL, no params object array"); // array ok
353  } else AliInfo("Do NOT recalibrate EMCAL, no params for pass"); // array pass ok
354  } else AliInfo("Do NOT recalibrate EMCAL, no params for run"); // run number array ok
355 
356  delete contRF;
357  } // Recalibration on
358 
359  // Energy Recalibration, apply on top of previous calibration factors
360  if(fRecoUtils->IsRunDepRecalibrationOn())
361  {
362  AliOADBContainer *contRFTD=new AliOADBContainer("");
363 
364  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
365 
366  TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
367 
368  //If it did not exist for this run, get closes one
369  if (!htd)
370  {
371  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
372  // let's get the closest runnumber instead then..
373  Int_t lower = 0;
374  Int_t ic = 0;
375  Int_t maxEntry = contRFTD->GetNumberOfEntries();
376 
377  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
378  lower = ic;
379  ic++;
380  }
381 
382  Int_t closest = lower;
383  if ( (ic<maxEntry) &&
384  (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
385  closest = ic;
386  }
387 
388  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
389  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
390  }
391 
392  if(htd)
393  {
394  AliInfo("Recalibrate (Temperature) EMCAL");
395 
396  for (Int_t ism=0; ism<nSM; ++ism)
397  {
398  for (Int_t icol=0; icol<48; ++icol)
399  {
400  for (Int_t irow=0; irow<24; ++irow)
401  {
402  Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
403 
404  Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
405  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
406  //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
407  // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
408  fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
409  } // columns
410  } // rows
411  } // SM loop
412  } else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
413 
414  delete contRFTD;
415  } // Run by Run T calibration
416 
417  // Time Recalibration
418  if(fRecoUtils->IsTimeRecalibrationOn())
419  {
420  AliOADBContainer *contTRF=new AliOADBContainer("");
421 
422  contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
423 
424  TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
425 
426  if(trecal)
427  {
428  // pass number should be pass1 except on Run1 and special cases
429  TString passM = pass;
430  if ( pass=="spc_calo" ) passM = "pass3";
431  if ( runnumber > 209121 ) passM = "pass1"; // run2 periods
432  if ( pass == "muon_calo_pass1" && runnumber > 209121 && runnumber < 244284 )
433  passM = "pass0";//period LHC15a-m
434 
435  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
436 
437  if(trecalpass)
438  {
439  AliInfo("Time Recalibrate EMCAL");
440  for (Int_t ibc = 0; ibc < 4; ++ibc)
441  {
442  TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
443 
444  if (h)
445  delete h;
446 
447  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
448 
449  if (!h)
450  {
451  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
452  continue;
453  }
454 
455  h->SetDirectory(0);
456 
457  fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
458  } // bunch crossing loop
459  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
460  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
461 
462  delete contTRF;
463  } // Time recalibration on
464 
465  // L1 Phase Time Recalibration
466  if( fRecoUtils->IsL1PhaseInTimeRecalibrationOn() )
467  {
468  // init default maps first
469  if (!fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationArray())
470  fRecoUtils->InitEMCALL1PhaseInTimeRecalibration() ;
471 
472  AliOADBContainer *contBC = new AliOADBContainer("");
473 
474  TFile *timeFile=new TFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"read");
475  if (!timeFile || timeFile->IsZombie())
476  {
477  AliFatal(Form("EMCALTimeL1PhaseCalib.root was not found in the path provided: %s",fOADBFilePath.Data()));
478  return ;
479  }
480 
481  if (timeFile) delete timeFile;
482 
483  contBC->InitFromFile(Form("%s/EMCALTimeL1PhaseCalib.root",fOADBFilePath.Data()),"AliEMCALTimeL1PhaseCalib");
484 
485  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
486  if (!arrayBC)
487  {
488  AliError(Form("No external L1 phase in time calibration set for run number: %d", runnumber));
489  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
490  }
491  else
492  {
493  // Only 1 L1 phase correction possible, except special cases
494  TString pass2 = "pass1";
495 
496  if ( pass=="muon_calo_pass1" && runnumber > 209121 && runnumber < 244284 )
497  pass2 = "pass0"; // period LHC15a-m
498 
499  TObjArray *arrayBCpass=(TObjArray*)arrayBC->FindObject(pass2);
500  if (!arrayBCpass)
501  {
502  AliError(Form("No external L1 phase in time calibration set for: %d -%s", runnumber,pass2.Data()));
503  fRecoUtils->SwitchOffL1PhaseInTimeRecalibration();
504  }
505  else AliInfo("Recalibrate L1 Phase time");
506 
507  if(arrayBCpass)
508  {
509  if ( DebugLevel()>0 ) arrayBCpass->Print();
510 
511  TH1C *h = fRecoUtils->GetEMCALL1PhaseInTimeRecalibrationForAllSM();
512  if (h) delete h;
513 
514  h = (TH1C*)arrayBCpass->FindObject(Form("h%d",runnumber));
515 
516  if (!h)
517  {
518  AliFatal(Form("There is no calibration histogram h%d for this run",runnumber));
519  return;
520  }
521 
522  h->SetDirectory(0);
523  fRecoUtils->SetEMCALL1PhaseInTimeRecalibrationForAllSM(h);
524  }
525  }
526 
527  delete contBC;
528  } // L1 Phase Time Recalibration
529 
530  // Parameters already set once, so do not it again
531  fOADBSet = kTRUE;
532 }
533 
534 //_________________________________________________
537 //_________________________________________________
539 {
540  fEvent = InputEvent();
541  if (!fEvent)
542  {
543  Warning("AccessOCDB","Event not available!!!");
544  return kFALSE;
545  }
546 
547  if (fEvent->GetRunNumber()==fRun)
548  return kTRUE;
549  fRun = fEvent->GetRunNumber();
550 
551  AliDebug(1,"Begin");
552 
553  AliCDBManager *cdb = AliCDBManager::Instance();
554 
555 
556  if (fOCDBpath.Length())
557  {
558  cdb->SetDefaultStorage(fOCDBpath.Data());
559  AliInfo(Form("Default storage %s",fOCDBpath.Data()));
560  }
561 
562  cdb->SetRun(fEvent->GetRunNumber());
563 
564  //
565  // EMCAL from RAW OCDB
566  if (fOCDBpath.Contains("alien:"))
567  {
568  cdb->SetSpecificStorage("EMCAL/Calib/Data","raw://");
569  cdb->SetSpecificStorage("EMCAL/Calib/Time","raw://");
570  cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","raw://");
571  }
572 
573  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
574 
575  return kTRUE;
576 }
577 
578 //_____________________________________________________
582 //_____________________________________________________
584 {
585  // Get the number of stored digits, to assign the index of the new ones
586  Int_t idigit = fDigitsArr->GetEntriesFast();
587 
588  Double_t fixTime = 615.*1e-9;
589 
590  for(Int_t j = 0; j < fgkNEMCalCells; j++)
591  {
592  // Newly created?
593  if ( !fTCardCorrCellsNew[j] ) continue;
594 
595  // Accept only if at least 10 MeV
596  if ( fTCardCorrCellsEner[j] < 0.01 ) continue;
597 
598  // Check if it was not masked
599  Float_t amp = 0;
600  Double_t time = 0;
601  if ( !fRecoUtils->AcceptCalibrateCell(j,0,amp,time,fCaloCells) ) continue;
602 
603 // printf("add new digit absId %d, accept? %d, digit %d, induced amp %2.2f\n",
604 // j,fRecoUtils->AcceptCalibrateCell(j,0,amp,time,fCaloCells),idigit,fTCardCorrCellsEner[j]);
605 
606  // Now add the cell to the digits list
607  //AliEMCALDigit* digit =
608  new((*fDigitsArr)[idigit]) AliEMCALDigit( -2, -2, j, fTCardCorrCellsEner[j], fixTime,AliEMCALDigit::kHG,idigit, 0, 0, 0);
609 
610  //digit->SetListOfParents(0,0x0,0x0); // not needed
611 
612  idigit++;
613  }// loop on all possible cells
614 }
615 
616 
617 //_____________________________________________________
622 //_____________________________________________________
624 {
625  fEvent = 0x0;
626 
627  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
628  Int_t eventN = Entry();
629  if(aodIH) eventN = aodIH->GetReadEntry();
630 
631  if (eventN > fMaxEvent)
632  return ;
633 
634  //printf("Clusterizer --- Event %d-- \n",eventN);
635 
636  // Check if input event are embedded events
637  // If so, take output event
638  if (aodIH && aodIH->GetMergeEvents())
639  {
640  fEvent = AODEvent();
641 
642  if(!aodIH->GetMergeEMCALCells())
643  AliFatal("Events merged but not EMCAL cells, check analysis settings!");
644 
645  AliDebug(1,"Use embedded events");
646 
647  AliDebug(1,Form("\t InputEvent N Clusters %d, N Cells %d",
648  InputEvent()->GetNumberOfCaloClusters(),InputEvent()->GetEMCALCells()->GetNumberOfCells()));
649 
650  AliDebug(1,Form("\t MergedEvent N Clusters %d, N Cells %d",
651  aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()));
652 
653 // if(DebugLevel() > 1)
654 // {
655 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
656 // {
657 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
658 // if(sigCluster->IsEMCAL()) AliInfo(Form("\t \t Signal cluster: i %d, E %f",icl,sigCluster->E()));
659 // }
660 // }
661 
662  AliDebug(1,Form("\t OutputEvent N Clusters %d, N Cells %d",
663  AODEvent()->GetNumberOfCaloClusters(), AODEvent()->GetEMCALCells()->GetNumberOfCells()));
664  }
665  else if(fInputFromFilter)
666  {
667  //printf("Get Input From Filtered AOD\n");
668  fEvent = AODEvent();
669  }
670  else
671  {
672  fEvent = InputEvent();
675  }
676 
677  if (!fEvent)
678  {
679  AliError("Event not available");
680  return ;
681  }
682 
683  //Recover the pointer to CaloCells container
684  fCaloCells = fEvent->GetEMCALCells();
685 
686  //Process events if there is a high energy cluster
687  if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
688 
689  //-------------------------------------------------------------------------------------
690  // Reject events if LED was firing, use only for LHC11a data
691  // Reject event if triggered by exotic cell and remove exotic cells if not triggered
692  //-------------------------------------------------------------------------------------
693 
694  if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
695 
696  if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
697 
698  //-------------------------------------------------------------------------------------
699  // Set the cluster array in the event (output or input)
700  //-------------------------------------------------------------------------------------
701 
702  if ( fFillAODFile )
703  {
704  //Magic line to write events to AOD filem put after event rejection
705  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
706  }
707  else if( !fOutputAODBranchSet )
708  {
709  // Create array and put it in the input event, if output AOD not selected, only once
710  InputEvent()->AddObject(fOutputAODBranch);
711  fOutputAODBranchSet = kTRUE;
712  AliInfo(Form("Add AOD branch <%s> to input event",fOutputAODBranchName.Data()));
713  }
714 }
715 
716 //____________________________________________________
721 //____________________________________________________
723 {
724  Int_t nClusters = fEvent->GetNumberOfCaloClusters();
725  Int_t nClustersOrg = 0;
726 
727  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
728  if(aodIH && aodIH->GetEventToMerge()) //Embedding
729  nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
730 
731  ResetArrays();
732 
733  // Loop on original clusters, get MC labels, cluster time (OLD AODs),
734  // or track matching residuals (if matching is not requested)
736  {
737  for (Int_t i = 0; i < nClusters; i++)
738  {
739  AliVCluster *clus = 0;
740  if(aodIH && aodIH->GetEventToMerge()) //Embedding
741  clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
742  else
743  clus = fEvent->GetCaloCluster(i);
744 
745  if(!clus) return;
746 
747  nClustersOrg++;
748 
749  if(!clus->IsEMCAL()) continue;
750 
751  Int_t label = clus->GetLabel();
752  Int_t label2 = -1 ;
753  if (clus->GetNLabels() >=2 ) label2 = clus->GetLabelAt(1) ;
754 
755  //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",
756  // i, clus->GetID(), clus->E(), clus->GetTOF()*1e9, clus->GetNCells(), clus->GetNLabels(), label);
757  //
758  //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++)
759  // printf("%d) Label %d, E dep frac %0.2f; ",
760  // imc, clus->GetLabelAt(imc),clus->GetClusterMCEdepFraction(imc));
761  //if(clus->GetNLabels() > 0) printf("\n");
762 
763  UShort_t * index = clus->GetCellsAbsId() ;
764  for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
765  {
766  fOrgClusterCellId[index[icell]] = i;
767  fCellTime[index[icell]] = clus->GetTOF();
768  fCellMatchdEta[index[icell]] = clus->GetTrackDz();
769  fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
770 
772  {
773  fCellLabels[index[icell]] = label;
774  fCellSecondLabels[index[icell]] = label2;
775  }
776  }
777  }
778  }
779 
780  // Do here induced cell energy assignation by T-Card correlation emulation, ONLY MC
782 
783  //
784  // Transform CaloCells into Digits
785  //
786  Int_t idigit = 0;
787  Int_t id = -1;
788  Float_t amp = -1;
789  Double_t time = -1;
790 
791  Int_t bc = InputEvent()->GetBunchCrossNumber();
792 
793  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
794  {
795  // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
796  id = fCaloCells->GetCellNumber(icell);
797  Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,fCaloCells);
798  time-=fConstantTimeShift*1e-9; // only in case of simulations done before 2015
799 
800  // Do not include cells with too low energy, nor exotic cell
801  if( amp < fRecParam->GetMinECut() ||
802  time > fRecParam->GetTimeMax() ||
803  time < fRecParam->GetTimeMin() ) accept = kFALSE;
804 
805  // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
807  {
808  time = fCellTime[id];
809  //printf("cell %d time org %f - ",id, time*1.e9);
810  fRecoUtils->RecalibrateCellTime(id,bc,time);
811  //printf("recal %f\n",time*1.e9);
812  }
813 
814  //Exotic?
815  if (accept && fRecoUtils->IsExoticCell(id,fCaloCells,bc))
816  accept = kFALSE;
817 
818  if( !accept )
819  {
820  AliDebug(2,Form("Remove channel absId %d, index %d of %d, amp %f, time %f",
821  id,icell, fCaloCells->GetNumberOfCells(), amp, time*1.e9));
822  continue;
823  }
824 
825  //
826  // MC
827  //
828 
829  // Old way
830  Int_t mcLabel = fCaloCells->GetMCLabel(icell);
831  Float_t eDep = amp;
832 
833  //printf("--- Selected cell %d--- amp %f\n",id,amp);
834 
835  // New way
836  TArrayI labeArr(0);
837  TArrayF eDepArr(0);
838  Int_t nLabels = 0;
839 
841  {
842  // Old way to recover/set the cell MC label
843  // Only possibility for old Run1 productions
844 
845  if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
846 
847  else if( fSetCellMCLabelFromCluster == 0 &&
849 
850  else mcLabel = -1; // found later
851 
852  // Last parameter of the digit object should be MC deposited energy,
853  // since it is not available in aliroot before year 2016, add just the cell amplitude so that
854  // we give more weight to the MC label of the cell with highest energy in the cluster
855 
856  Float_t efrac = fCaloCells->GetEFraction(icell);
857 
858  // When checking the MC of digits, give weight to cells with embedded signal
859  if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
860 
861  eDep *= efrac ;
862  }
863  else // fSetCellMCLabelFromEdepFrac = true
864  {
865  // New way, valid only for MC productions with aliroot > v5-07-21
866  mcLabel = -1;
867 
868  // Map the digit to cell index for later to calculate the cell MC energy deposition map
869  fCellLabels[id] = idigit;
870  //printf("\t absId %d, idigit %d\n",id,idigit);
871 
872  if(fOrgClusterCellId[id] < 0) continue; // index can be negative if noisy cell that did not form cluster
873 
874  AliVCluster *clus = 0;
875  Int_t iclus = fOrgClusterCellId[id];
876 
877  if(iclus < 0)
878  {
879  AliInfo("Negative original cluster index, skip \n");
880  continue;
881  }
882 
883  if(aodIH && aodIH->GetEventToMerge()) //Embedding
884  clus = aodIH->GetEventToMerge()->GetCaloCluster(iclus); //Get clusters directly from embedded signal
885  else
886  clus = fEvent->GetCaloCluster(iclus);
887 
888  fRecoUtils->RecalculateCellLabelsRemoveAddedGenerator(id, clus, MCEvent(), amp, labeArr, eDepArr);
889  nLabels = labeArr.GetSize();
890 
891  } // cell MC label, new
892 
893  // Apply here the found induced energies
894  if ( fTCardCorrEmulation )
895  {
896 // if( TMath::Abs(fTCardCorrCellsEner[id]) > 0.001 )
897 // printf("add energy to digit %d, absId %d: amp %2.2f + %2.2f\n",idigit,id,amp,fTCardCorrCellsEner[id]);
898  amp+=fTCardCorrCellsEner[id];
899  }
900 
901  //
902  // Create the digit
903  //
904  if(amp <= 0.01) continue ; // accept if > 10 MeV
905 
906  //if(amp > 1.5)printf("*** Add digit *** amp %f, nlabels %d, label %d, found %d, edep fract tot %f, ncluslabels %d\n",amp, nLabels,mcLabel,found,edepTotFrac,nClusLabels);
907 
908  AliEMCALDigit* digit = new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, eDep);
909 
910  if(nLabels > 0)
911  digit->SetListOfParents(nLabels,labeArr.GetArray(),eDepArr.GetArray());
912 
913  idigit++;
914 
915  }
916 
917  fDigitsArr->Sort();
918 
919  // Call after Sort, if not it screws up digits index order in clusterization
921 
922 
923  //-------------------------------------------------------------------------------------
924  // Do the clusterization
925  //-------------------------------------------------------------------------------------
926 
927  fClusterizer->Digits2Clusters("");
928 
929  //-------------------------------------------------------------------------------------
930  // Transform the recpoints into AliVClusters
931  //-------------------------------------------------------------------------------------
932 
934 
935  if(!fCaloClusterArr)
936  {
937  AliWarning(Form("No array with CaloClusters, input RecPoints entries %d",fClusterArr->GetEntriesFast()));
938  return;
939  }
940 
941  AliDebug(1,Form("N clusters: before recluster %d, after recluster %d, recpoints %d",
942  nClustersOrg, fCaloClusterArr->GetEntriesFast(),fClusterArr->GetEntriesFast()));
943 
944 // if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
945 // {
946 // AliInfo("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
947 // fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
948 // fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
949 // }
950 }
951 
952 
953 //_____________________________________________________
955 //_____________________________________________________
957 {
958  Double_t cellAmplitude = 0;
959  Double_t cellTime = 0;
960  Short_t cellNumber = 0;
961  Int_t cellMCLabel = 0;
962  Double_t cellEFrac = 0;
963  Int_t nClustersOrg = 0;
964 
965  // Fill the array with the EMCAL clusters, copy them
966  for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
967  {
968  AliVCluster *clus = fEvent->GetCaloCluster(i);
969  if(clus->IsEMCAL())
970  {
971  //recalibrate/remove bad channels/etc if requested
972  if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
973  {
974  continue;
975  }
976 
977  if(fRecoUtils->IsRecalibrationOn())
978  {
979  //Calibrate cluster
980  fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, fCaloCells);
981 
982  //CalibrateCells
983  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
984  {
985  if (fCaloCells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
986  break;
987 
988  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
989  fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
990  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
991 
992  //Do not include bad channels found in analysis?
993  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
994  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
995  continue;
996 
997  fCaloCells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
998 
999  }// cells loop
1000  }// recalibrate
1001 
1002  //Cast to ESD or AOD, needed to create the cluster array
1003  AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
1004  AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
1005 
1006  if (esdCluster)
1007  {
1008  fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
1009  }//ESD
1010  else if(aodCluster)
1011  {
1012  fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
1013  }//AOD
1014  else
1015  AliWarning("Wrong CaloCluster type?");
1016 
1017  nClustersOrg++;
1018  }
1019  }
1020 
1021  //Do the unfolding
1022  fUnfolder->UnfoldClusters(fCaloClusterArr, fCaloCells);
1023 
1024  //CLEAN-UP
1025  fUnfolder->Clear();
1026 }
1027 
1028 //_____________________________________________________
1030 //_____________________________________________________
1032 {
1033  AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
1034  Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
1035 
1036  AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
1037  aodEMcells.CreateContainer(nEMcell);
1038  aodEMcells.SetType(AliVCaloCells::kEMCALCell);
1039  Double_t calibFactor = 1.;
1040  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
1041  {
1042  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1043  fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
1044  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1045 
1046  if(fRecoUtils->IsRecalibrationOn())
1047  {
1048  calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
1049  }
1050 
1051  if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
1052  { //Channel is not declared as bad
1053  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
1054  eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
1055  }
1056  else
1057  {
1058  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
1059  }
1060  }
1061 
1062  aodEMcells.Sort();
1063 }
1064 
1065 //__________________________________________________
1067 //__________________________________________________
1069 {
1070  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
1071  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
1072 
1073  Double_t pos[3] ;
1074  Double_t covVtx[6];
1075  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
1076 
1077  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
1078  if(!header) AliFatal("Not a standard AOD");
1079  header->SetRunNumber(fEvent->GetRunNumber());
1080 
1081  if(esdevent)
1082  {
1083  TTree* tree = fInputHandler->GetTree();
1084  if (tree)
1085  {
1086  TFile* file = tree->GetCurrentFile();
1087  if (file) header->SetESDFileName(file->GetName());
1088  }
1089  }
1090  else if (aodevent) {
1091  AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
1092  if(!aodheader) AliFatal("Not a standard AOD");
1093  header->SetESDFileName(aodheader->GetESDFileName());
1094  }
1095 
1096  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
1097  header->SetOrbitNumber(fEvent->GetOrbitNumber());
1098  header->SetPeriodNumber(fEvent->GetPeriodNumber());
1099  header->SetEventType(fEvent->GetEventType());
1100 
1101  // Centrality
1102  if(fUseAliCentrality)
1103  {
1104  if(fEvent->GetCentrality())
1105  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
1106  else
1107  header->SetCentrality(0);
1108  }
1109 
1110  // Trigger
1111  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
1112 
1113  header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
1114 
1115  header->SetTriggerMask(fEvent->GetTriggerMask());
1116  header->SetTriggerCluster(fEvent->GetTriggerCluster());
1117 
1118  if (esdevent)
1119  {
1120  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
1121  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
1122  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
1123  }
1124  else if (aodevent)
1125  {
1126  header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
1127  header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
1128  header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
1129  }
1130 
1131  header->SetMagneticField(fEvent->GetMagneticField());
1132  //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
1133 
1134  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
1135  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
1136  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
1137  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
1138  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
1139 
1140  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
1141  Float_t diamcov[3];
1142  fEvent->GetDiamondCovXY(diamcov);
1143  header->SetDiamond(diamxy,diamcov);
1144  if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
1145  else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
1146 
1147  //
1148  Int_t nVertices = 1 ;/* = prim. vtx*/;
1149  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
1150 
1151  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
1152 
1153  // Access to the AOD container of vertices
1154  TClonesArray &vertices = *(AODEvent()->GetVertices());
1155  Int_t jVertices=0;
1156 
1157  // Add primary vertex. The primary tracks will be defined
1158  // after the loops on the composite objects (V0, cascades, kinks)
1159  fEvent->GetPrimaryVertex()->GetXYZ(pos);
1160  Float_t chi = 0;
1161  if (esdevent)
1162  {
1163  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1164  chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
1165  }
1166  else if (aodevent)
1167  {
1168  aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
1169  chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
1170  }
1171 
1172  AliAODVertex * primary = new(vertices[jVertices++])
1173  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
1174  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
1175  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
1176 }
1177 
1178 //___________________________________________________________
1182 //___________________________________________________________
1184 {
1185  // First recalculate track-matching for the new clusters
1186  if(fDoTrackMatching)
1187  {
1188  fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom,MCEvent());
1189  }
1190 
1191  // Put the new clusters in the AOD list
1192 
1193  Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
1194 
1195  for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
1196  {
1197  AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
1198 
1199  newCluster->SetID(i);
1200 
1201  // Correct cluster energy non linearity
1202  newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
1203 
1204  //Add matched track
1205  if(fDoTrackMatching)
1206  {
1207  Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
1208  if(trackIndex >= 0)
1209  {
1210  newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
1211  AliDebug(2,Form("Matched Track index %d to new cluster %d",trackIndex,i));
1212  }
1213 
1214  Float_t dR = 999., dZ = 999.;
1215  fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dZ,dR);
1216  newCluster->SetTrackDistance(dR,dZ);
1217  }
1218  else
1219  {// Assign previously assigned matched track in reco, very very rough
1220  Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
1221  newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
1222  }
1223 
1224  //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
1225  //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
1226  //printf("\n");
1227  //
1228  //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",
1229  // i, newCluster->GetID(), newCluster->E(), newCluster->GetTOF()*1e9, newCluster->GetNCells(), newCluster->GetNLabels(), newCluster->GetLabel());
1230  //
1231  //for(Int_t imc = 0; imc < newCluster->GetNLabels(); imc++)
1232  // printf("%d) Label %d, E dep frac %0.2f; ",
1233  // imc,newCluster->GetLabelAt(imc),newCluster->GetClusterMCEdepFraction(imc));
1234  //if(newCluster->GetNLabels() > 0) printf("\n");
1235 
1236  // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
1237  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fCaloCells, newCluster);
1238 
1239  new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
1240 
1241  AliDebug(2,Form("New cluster %d of %d, energy %f, mc label %d",
1242  newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()));
1243 
1244  } // cluster loop
1245 
1246  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
1247 }
1248 
1249 
1250 //________________________________________________________________
1252 //________________________________________________________________
1254 {
1255  if(fUseAliCentrality)
1256  {
1257  if(GetCentrality()) return GetCentrality()->GetCentralityPercentile(fCentralityClass) ;
1258  else return -1. ;
1259  }
1260  else
1261  {
1262  if(GetMultSelCen()) return GetMultSelCen()->GetMultiplicityPercentile(fCentralityClass,kTRUE) ;
1263  else return -1. ;
1264  }
1265 
1266  return -1.;
1267 }
1268 
1269 //_______________________________________________
1271 //_______________________________________________
1273 {
1274  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
1275  {
1276  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null");
1277  return TString("");
1278  }
1279 
1280  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
1281  {
1282  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null");
1283  return TString("");
1284  }
1285 
1286  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1287  if (pass.Contains("ass1")) return TString("pass1");
1288  else if (pass.Contains("ass2")) return TString("pass2");
1289  else if (pass.Contains("ass3")) return TString("pass3");
1290  else if (pass.Contains("ass4")) return TString("pass4");
1291  else if (pass.Contains("ass5")) return TString("pass5");
1292  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1293  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1294  {
1295  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1296  return TString("pass1");
1297  }
1298  else if (pass.Contains("LHC14a1a"))
1299  {
1300  AliInfo("Enable EMCal energy calibration for this MC production!!");
1301 
1302  fRecoUtils->SwitchOnRecalibration();
1303 
1304  return TString("LHC14a1a");
1305  }
1306 
1307  // No condition fullfilled, give a default value
1308  AliInfo("Pass number string not found");
1309 
1310  return TString("");
1311 }
1312 
1313 //_________________________________________
1316 //_________________________________________
1318 {
1319  if(fDebug >=0) (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1320 
1321  fOADBSet = kFALSE;
1322  if(fOADBFilePath == "") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
1323 
1324  fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1325 
1326  if(!fRecParam) fRecParam = new AliEMCALRecParam;
1327  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1328 
1329  if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1330  if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1331  if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1332  fRejectBelowThreshold = kFALSE;
1333 
1334  // Centrality
1335  if(fCentralityClass == "") fCentralityClass = "V0M";
1336 
1337  if (fOCDBpath == "") fOCDBpath = "raw://" ;
1338  if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1339 
1340  if(gROOT->LoadMacro(fConfigName) >=0)
1341  {
1342  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
1343  AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1344  fGeomName = clus->fGeomName;
1346  fOCDBpath = clus->fOCDBpath;
1347  fAccessOCDB = clus->fAccessOCDB;
1348  fRecParam = clus->fRecParam;
1349  fJustUnfold = clus->fJustUnfold;
1350  fFillAODFile = clus->fFillAODFile;
1351  fRecoUtils = clus->fRecoUtils;
1352  fConfigName = clus->fConfigName;
1353  fMaxEvent = clus->fMaxEvent;
1355  fUpdateCell = clus->fUpdateCell;
1357  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1359  fCentralityBin[0] = clus->fCentralityBin[0];
1360  fCentralityBin[1] = clus->fCentralityBin[1];
1362  }
1363 }
1364 
1365 //_______________________________________________________
1368 //_______________________________________________________
1370 {
1371  if (fJustUnfold)
1372  {
1373  // init the unfolding afterburner
1374  delete fUnfolder;
1375  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1376  return;
1377  }
1378 
1379  //First init the clusterizer
1380  delete fClusterizer;
1381  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1382  fClusterizer = new AliEMCALClusterizerv1 (fGeom);
1383  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1384  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1385  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1386  {
1387  fClusterizer = new AliEMCALClusterizerNxN(fGeom);
1388  fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1389  fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1390  }
1391  else
1392  {
1393  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1394  }
1395 
1396  //Now set the parameters
1397  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1398  fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1399  fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1400  fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1401  fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1402  fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1403  fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1404  fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1405  fClusterizer->SetTimeCalibration ( fRecParam->IsTimeCalibrationOn() );
1406  fClusterizer->SetInputCalibrated ( kTRUE );
1407  fClusterizer->SetJustClusters ( kTRUE );
1408 
1409  // Initialize the cluster rec points and digits arrays and get them.
1410  fClusterizer->SetOutput(0);
1411  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1412  fDigitsArr = fClusterizer->GetDigits();
1413 
1414  //In case of unfolding after clusterization is requested, set the corresponding parameters
1415  if(fRecParam->GetUnfold())
1416  {
1417  Int_t i=0;
1418  for (i = 0; i < 8; i++)
1419  {
1420  fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1421  }//end of loop over parameters
1422 
1423  for (i = 0; i < 3; i++)
1424  {
1425  fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1426  fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1427  }//end of loop over parameters
1428  fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1429  fClusterizer->InitClusterUnfolding();
1430 
1431  }// to unfold
1432 }
1433 
1434 //_________________________________________________
1438 //_________________________________________________
1440 {
1441  if(fGeomMatrixSet) return;
1442 
1443  Int_t runnumber = InputEvent()->GetRunNumber() ;
1444  if (!fGeom)
1445  {
1446  if(fGeomName=="")
1447  {
1448  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(runnumber);
1449  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fGeom->GetName(),runnumber));
1450  }
1451  else
1452  {
1453  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1454  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeomName.Data()));
1455  }
1456 
1457  // Init geometry, I do not like much to do it like this ...
1458  if(fImportGeometryFromFile && !gGeoManager)
1459  {
1460  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1461  {
1462  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1463  if (runnumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1464  else if (runnumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1465  else if (runnumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1466  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >=2015
1467  }
1468 
1469  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1470 
1471  TGeoManager::Import(fImportGeometryFilePath) ;
1472  }
1473 
1474  AliDebug(1,Form("Init for run=%d",runnumber));
1475  if (!gGeoManager) AliDebug(1,"Careful!, gGeoManager not loaded, load misalign matrices");
1476  } // geometry pointer did not exist before
1477 
1478  if(fLoadGeomMatrices)
1479  {
1480  AliInfo("Load user defined EMCAL geometry matrices");
1481 
1482  // OADB if available
1483  AliOADBContainer emcGeoMat("AliEMCALgeo");
1484  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1485  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1486 
1487  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1488  {
1489 
1490  if (!fGeomMatrix[mod]) // Get it from OADB
1491  {
1492  AliDebug(2,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
1493  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1494 
1495  fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1496  }
1497 
1498  if(fGeomMatrix[mod])
1499  {
1500  if(DebugLevel() > 1)
1501  fGeomMatrix[mod]->Print();
1502 
1503  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1504  }
1505  else if(gGeoManager)
1506  {
1507  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
1508  fGeom->SetMisalMatrix(fGeom->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
1509  }
1510  else
1511  {
1512  AliError(Form("Alignment atrix for SM %d is not available",mod));
1513  }
1514 
1515  fGeomMatrixSet=kTRUE;
1516 
1517  }//SM loop
1518  }//Load matrices
1519  else if(!gGeoManager)
1520  {
1521  AliInfo("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1522  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1523  if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1524  {
1525  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
1526  }//AOD
1527  else
1528  {
1529  AliDebug(1,"Load Misaligned matrices");
1530 
1531  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1532 
1533  if(!esd)
1534  {
1535  AliError("This event does not contain ESDs?");
1536  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1537  return;
1538  }
1539 
1540  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1541  {
1542  if(DebugLevel() > 1)
1543  esd->GetEMCALMatrix(mod)->Print();
1544 
1545  if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1546 
1547  }
1548 
1549  fGeomMatrixSet=kTRUE;
1550 
1551  } // ESD
1552  } // Load matrices from Data
1553 }
1554 
1555 //____________________________________________________
1560 //____________________________________________________
1562 {
1563  if(!fRemoveExoticEvents) return kFALSE;
1564 
1565  // Loop on cells
1566 
1567  Float_t totCellE = 0;
1568  Int_t bc = InputEvent()->GetBunchCrossNumber();
1569 
1570  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1571  {
1572  Float_t ecell = 0 ;
1573  Double_t tcell = 0 ;
1574 
1575  Int_t absID = fCaloCells->GetCellNumber(icell);
1576  Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,fCaloCells);
1577  tcell-=fConstantTimeShift*1e-9;// Only for MC simulations done before 2015
1578 
1579  if(accept && !fRecoUtils->IsExoticCell(absID,fCaloCells,bc)) totCellE += ecell;
1580  }
1581 
1582  // TString triggerclasses = event->GetFiredTriggerClasses();
1583  // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1584  // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1585  // return;
1586  //
1587 
1588  //printf("TotE cell %f\n",totCellE);
1589  if(totCellE < 1) return kTRUE;
1590 
1591  return kFALSE;
1592 }
1593 
1594 //________________________________________________________________
1596 //________________________________________________________________
1598 {
1599  if(!fRemoveLEDEvents) return kFALSE;
1600 
1601  //check events of LHC11a period
1602  if(run < 146858 || run > 146860) return kFALSE ;
1603 
1604  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1605  Int_t ncellsSM3 = 0;
1606  for(Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1607  {
1608  if ( fCaloCells->GetAmplitude (icell) > 0.1 &&
1609  fCaloCells->GetCellNumber(icell)/(24*48)==3 ) ncellsSM3++;
1610  }
1611 
1612  TString triggerclasses = fEvent->GetFiredTriggerClasses();
1613 
1614  Int_t ncellcut = 21;
1615  if(triggerclasses.Contains("EMC")) ncellcut = 35;
1616 
1617  if( ncellsSM3 >= ncellcut)
1618  {
1619  AliInfo(Form("Reject event %d with ncells in SM3 %d",(Int_t)Entry(),ncellsSM3));
1620  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1621  return kTRUE;
1622  }
1623 
1624  return kFALSE;
1625 }
1626 
1627 
1628 //_______________________________________________________
1631 //_______________________________________________________
1633 {
1634  Int_t id = -1;
1635  Float_t amp = -1;
1636 
1637  // Loop on all cells with signal
1638  for (Int_t icell = 0; icell < fCaloCells->GetNumberOfCells(); icell++)
1639  {
1640  id = fCaloCells->GetCellNumber(icell);
1641  amp = fCaloCells->GetAmplitude (icell); // fCaloCells->GetCellAmplitude(id);
1642 
1643  //
1644  // First get the SM, col-row of this tower
1645  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
1646  fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
1647  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
1648 
1649  //
1650  // Determine randomly if we want to create a correlation for this cell,
1651  // depending the SM number of the cell
1652  Float_t rand = fRandom.Uniform(0, 1);
1653 
1654 // printf("icell %d, absId %d, ism %d, icol %d, irow %d, random %2.2f > %2.2f ?; \n \t amp %2.3f added before %2.3f\n",
1655 // icell,id,imod,ieta,iphi,rand,fTCardCorrInduceEnerProb[imod],amp,fTCardCorrCellsEner[id]);
1656 
1657  if ( rand > fTCardCorrInduceEnerProb[imod] ) continue;
1658 
1659  //
1660  // Get the absId of the cells in the cross and same T-Card
1661  Int_t absIDup = -1;
1662  Int_t absIDdo = -1;
1663  Int_t absIDlr = -1;
1664  Int_t absIDuplr = -1;
1665  Int_t absIDdolr = -1;
1666 
1667  Int_t absIDup2 = -1;
1668  Int_t absIDup2lr = -1;
1669  Int_t absIDdo2 = -1;
1670  Int_t absIDdo2lr = -1;
1671 
1672  // Only 2 columns in the T-Card, +1 for even and -1 for odd with respect reference cell
1673  Int_t colShift = 0;
1674  if ( (ieta%2) && ieta <= AliEMCALGeoParams::fgkEMCALCols-1 ) colShift = -1;
1675  if ( !(ieta%2) && ieta >= 0 ) colShift = +1;
1676 
1677  absIDlr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+colShift);
1678 
1679  // Check up / down cells from reference cell not out of SM and in same T-Card
1680  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-1 )
1681  {
1682  absIDup = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta);
1683  absIDuplr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+1, ieta+colShift);
1684  }
1685 
1686  if ( iphi > 0 )
1687  {
1688  absIDdo = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta);
1689  absIDdolr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-1, ieta+colShift);
1690  }
1691 
1692  // Check 2 up / 2 down cells from reference cell not out of SM
1693  if ( iphi < AliEMCALGeoParams::fgkEMCALRows-2 )
1694  {
1695  absIDup2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta);
1696  absIDup2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi+2, ieta+colShift);
1697  }
1698 
1699  if ( iphi > 1 )
1700  {
1701  absIDdo2 = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta);
1702  absIDdo2lr = fGeom->GetAbsCellIdFromCellIndexes(imod, iphi-2, ieta+colShift);
1703  }
1704 
1705  // In same T-Card?
1706  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+1)/8) ) { absIDup = -1 ; absIDuplr = -1 ; }
1707  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-1)/8) ) { absIDdo = -1 ; absIDdolr = -1 ; }
1708  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi+2)/8) ) { absIDup2 = -1 ; absIDup2lr = -1 ; }
1709  if ( TMath::FloorNint(iphi/8) != TMath::FloorNint((iphi-2)/8) ) { absIDdo2 = -1 ; absIDdo2lr = -1 ; }
1710 
1711 
1712 // printf("\t cross absId: r1 %d, r2 %d, c %d\n",absIDup,absIDdo,absIDlr);
1713 // if ( absIDlr < 0 ) printf( "\t \t *** ism %d, iphi %d, ieta %d; absId(ieta+1) %d, absId(ieta-1) %d, mod eta %d, max eta %d\n",imod,iphi,ieta,
1714 // fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta+1),fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, ieta-1),
1715 // ieta%2, AliEMCALGeoParams::fgkEMCALCols-1);
1716 // if ( absIDup < 0 ) printf( "\t \t *** ism %d, iphi %d, ieta %d; absId(iphi+1) %d\n",imod,iphi,ieta,
1717 // fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, iphi+1));
1718 //
1719 // if ( absIDdo < 0 ) printf( "\t \t *** ism %d, iphi %d, ieta %d; absId(iphi-1) %d\n",imod,iphi,ieta,
1720 // fGeom->GetAbsCellIdFromCellIndexes(imod, iphi, iphi-1));
1721 
1722  //
1723  // Check if they are not declared bad or exist
1724  Bool_t okup = AcceptCell(absIDup );
1725  Bool_t okdo = AcceptCell(absIDdo );
1726  Bool_t oklr = AcceptCell(absIDlr );
1727  Bool_t okuplr = AcceptCell(absIDuplr );
1728  Bool_t okdolr = AcceptCell(absIDdolr );
1729  Bool_t okup2 = AcceptCell(absIDup2 );
1730  Bool_t okdo2 = AcceptCell(absIDdo2 );
1731  Bool_t okup2lr= AcceptCell(absIDup2lr);
1732  Bool_t okdo2lr= AcceptCell(absIDdo2lr);
1733 
1734 // if(amp > 2)
1735 // printf("ecell %2.2f, amp %2.2f: F1 %2.2e, F2 %2.2e, F3 %2.2e; w1 %2.2e w2 %2.2e w3 %2.2e \n",ecell,amp,
1736 // fTCardCorrInduceEnerFrac[0]+amp*fTCardCorrInduceEnerFracP1[0],
1737 // fTCardCorrInduceEnerFrac[1]+amp*fTCardCorrInduceEnerFracP1[1],
1738 // fTCardCorrInduceEnerFrac[1]+amp*fTCardCorrInduceEnerFracP1[2],
1739 // fTCardCorrInduceEnerFracWidth[0],fTCardCorrInduceEnerFracWidth[1],fTCardCorrInduceEnerFracWidth[2]);
1740  //
1741  // Generate some energy for the nearby cells in same TCard , depending on this cell energy
1742  // Check if originally the tower had no or little energy, in which case tag it as new
1744  Float_t fracupdownleri = fRandom.Gaus(fTCardCorrInduceEnerFrac[1]+amp*fTCardCorrInduceEnerFracP1[1],fTCardCorrInduceEnerFracWidth[1]);
1745  Float_t fracleri = fRandom.Gaus(fTCardCorrInduceEnerFrac[2]+amp*fTCardCorrInduceEnerFracP1[2],fTCardCorrInduceEnerFracWidth[2]);
1746  Float_t frac2nd = fRandom.Gaus(fTCardCorrInduceEnerFrac[3]+amp*fTCardCorrInduceEnerFracP1[3],fTCardCorrInduceEnerFracWidth[3]);
1747 
1748 // // CAREFUL: <<Apply the same added shift to all>>.
1749 // Float_t fracCupdown = fracRupdown;
1750 // Float_t fracC = fracRupdown;
1751 
1752 // printf("\t Added energy (fractions %2.3f-%2.3f-%2.3f) -> gaus (%2.3f-%2.3f-%2.3f) on top of previous %2.3f: \n",
1753 // fTCardCorrInduceEnerFrac[0], fTCardCorrInduceEnerFrac[1], fTCardCorrInduceEnerFrac[2],
1754 // fracRupdown,fracC,fracCupdown, fTCardCorrCellsEner[id]);
1755 
1756  if ( okup )
1757  {
1758  fTCardCorrCellsEner[absIDup] += amp*fracupdown;
1759  if ( fCaloCells->GetCellAmplitude(absIDup) < 0.01 ) fTCardCorrCellsNew[absIDup] = kTRUE;
1760 // printf("\t \t r1 absId %d induced amp %2.3f, accumulated %2.3f (new %d); \n",
1761 // absIDup,amp*fracRupdown,fTCardCorrCellsEner[absIDup],fTCardCorrCellsNew[absIDup]);
1762  }
1763 
1764  if ( okdo )
1765  {
1766  fTCardCorrCellsEner[absIDdo] += amp*fracupdown;
1767  if ( fCaloCells->GetCellAmplitude(absIDdo) < 0.01 ) fTCardCorrCellsNew[absIDdo] = kTRUE;
1768 // printf("\t \t r2 absId %d induced amp %2.3f, accumulated %2.3f (new %d);\n ",
1769 // absIDdo,amp*fracRupdown,fTCardCorrCellsEner[absIDdo],fTCardCorrCellsNew[absIDdo]);
1770  }
1771 
1772  if ( oklr )
1773  {
1774  fTCardCorrCellsEner[absIDlr] += amp*fracleri;
1775  if ( fCaloCells->GetCellAmplitude(absIDlr) < 0.01 ) fTCardCorrCellsNew[absIDlr] = kTRUE;
1776 // printf("\t \t c absId %d induced amp %2.3f, accumulated %2.3f (new %d). \n",
1777 // absIDlr,amp*fracC,fTCardCorrCellsEner[absIDlr],fTCardCorrCellsNew[absIDlr]);
1778  }
1779 
1780  if ( okuplr )
1781  {
1782  fTCardCorrCellsEner[absIDuplr] += amp*fracupdownleri;
1783  if ( fCaloCells->GetCellAmplitude(absIDuplr ) < 0.01 ) fTCardCorrCellsNew[absIDuplr] = kTRUE;
1784  // printf("\t \t c absId %d induced amp %2.3f, accumulated %2.3f (new %d). \n",
1785  // absIDuplr,amp*fracCupdown,fTCardCorrCellsEner[absIDuplr],fTCardCorrCellsNew[absIDuplr]);
1786  }
1787 
1788  if ( okdolr )
1789  {
1790  fTCardCorrCellsEner[absIDdolr] += amp*fracupdownleri;
1791  if ( fCaloCells->GetCellAmplitude(absIDdolr ) < 0.01 ) fTCardCorrCellsNew[absIDdolr] = kTRUE;
1792  // printf("\t \t c absId %d induced amp %2.3f, accumulated %2.3f (new %d). \n",
1793  // absIDdolr,amp*fracCupdown,fTCardCorrCellsEner[absIDdolr],fTCardCorrCellsNew[absIDdolr]);
1794  }
1795 
1796 
1797  if ( okup2 )
1798  {
1799  fTCardCorrCellsEner[absIDup2] += amp*frac2nd;
1800  if ( fCaloCells->GetCellAmplitude(absIDup2) < 0.01 ) fTCardCorrCellsNew[absIDup2] = kTRUE;
1801  // printf("\t \t r1 absId %d induced amp %2.3f, accumulated %2.3f (new %d); \n",
1802  // absIDup2,amp*fracRupdown,fTCardCorrCellsEner[absIDup2],fTCardCorrCellsNew[absIDup2]);
1803  }
1804 
1805  if ( okup2lr )
1806  {
1807  fTCardCorrCellsEner[absIDup2lr] += amp*frac2nd;
1808  if ( fCaloCells->GetCellAmplitude(absIDup2lr) < 0.01 ) fTCardCorrCellsNew[absIDup2lr] = kTRUE;
1809  // printf("\t \t r2 absId %d induced amp %2.3f, accumulated %2.3f (new %d);\n ",
1810  // absIDr2,amp*fracRupdown,fTCardCorrCellsEner[absIDr2],fTCardCorrCellsNew[absIDr2]);
1811  }
1812 
1813  if ( okdo2 )
1814  {
1815  fTCardCorrCellsEner[absIDdo2] += amp*frac2nd;
1816  if ( fCaloCells->GetCellAmplitude(absIDdo2) < 0.01 ) fTCardCorrCellsNew[absIDdo2] = kTRUE;
1817  // printf("\t \t r1 absId %d induced amp %2.3f, accumulated %2.3f (new %d); \n",
1818  // absIDdo2,amp*fracRupdown,fTCardCorrCellsEner[absIDdo2],fTCardCorrCellsNew[absIDdo2]);
1819  }
1820 
1821  if ( okdo2lr )
1822  {
1823  fTCardCorrCellsEner[absIDdo2lr] += amp*frac2nd;
1824  if ( fCaloCells->GetCellAmplitude(absIDdo2lr) < 0.01 ) fTCardCorrCellsNew[absIDdo2lr] = kTRUE;
1825  // printf("\t \t r2 absId %d induced amp %2.3f, accumulated %2.3f (new %d);\n ",
1826  // absIDr2,amp*fracRupdown,fTCardCorrCellsEner[absIDr2],fTCardCorrCellsNew[absIDr2]);
1827  }
1828 
1829 
1830  //
1831  // Subtract the added energy to main cell, if energy conservation is requested
1833  {
1834  if ( oklr ) fTCardCorrCellsEner[id] -= amp*fracleri;
1835  if ( okuplr ) fTCardCorrCellsEner[id] -= amp*fracupdownleri;
1836  if ( okdolr ) fTCardCorrCellsEner[id] -= amp*fracupdownleri;
1837  if ( okup ) fTCardCorrCellsEner[id] -= amp*fracupdown;
1838  if ( okdo ) fTCardCorrCellsEner[id] -= amp*fracupdown;
1839  if ( okup2 ) fTCardCorrCellsEner[id] -= amp*frac2nd;
1840  if ( okup2lr ) fTCardCorrCellsEner[id] -= amp*frac2nd;
1841  if ( okdo2 ) fTCardCorrCellsEner[id] -= amp*frac2nd;
1842  if ( okdo2lr ) fTCardCorrCellsEner[id] -= amp*frac2nd;
1843  //printf("\t conserve energy, remove %2.3f, from %d\n",fTCardCorrCellsEner[id],id);
1844  } // conserve energy
1845 
1846  } // cell loop
1847 
1848 }
1849 
1850 //_______________________________________________________
1854 //_______________________________________________________
1856 {
1857  Int_t j = 0;
1858  for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1859  {
1860  AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1861 
1862  const Int_t ncells = recPoint->GetMultiplicity();
1863  Int_t ncellsTrue = 0;
1864 
1865  if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1866 
1867  // cells and their amplitude fractions
1868  UShort_t absIds[ncells];
1869  Double32_t ratios[ncells];
1870 
1871  //For later check embedding
1872  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1873 
1874  Float_t clusterE = 0;
1875  for (Int_t c = 0; c < ncells; c++)
1876  {
1877  AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1878 
1879  absIds[ncellsTrue] = digit->GetId();
1880  ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1881 
1882  if ( !fRecParam->GetUnfold() && (ratios[ncellsTrue] > 1 || ratios[ncellsTrue] < 1) )
1883  AliWarning(Form("recpoint cell E %2.3f but digit E %2.3f and no unfolding", recPoint->GetEnergiesList()[c], digit->GetAmplitude()));
1884 
1885  // In case of unfolding, remove digits with unfolded energy too low
1886  if(fSelectCell)
1887  {
1888  if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1889  {
1890 
1891  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
1892  recPoint->GetEnergiesList()[c],digit->GetAmplitude()));
1893  continue;
1894  } // if cuts
1895  }// Select cells
1896 
1897  //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1898  clusterE +=recPoint->GetEnergiesList()[c];
1899 
1900  // In case of embedding, fill ratio with amount of signal,
1901  if (aodIH && aodIH->GetMergeEvents())
1902  {
1903  //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1904  AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1905  AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1906 
1907  Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1908  //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1909  Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1910  //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1911 
1912  if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1913  //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1914 
1915  }//Embedding
1916 
1917  ncellsTrue++;
1918 
1919  }// cluster cell loop
1920 
1921  if (ncellsTrue < 1)
1922  {
1923  AliDebug(2,Form("Skipping cluster with no cells avobe threshold E = %f, ncells %d",
1924  recPoint->GetEnergy(), ncells));
1925  continue;
1926  }
1927 
1928  //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1929 
1930  if(clusterE < fRecParam->GetClusteringThreshold())
1931  {
1932  AliDebug(2,Form("Remove cluster with energy below seed threshold %f",clusterE));
1933  continue;
1934  }
1935 
1936  TVector3 gpos;
1937  Float_t g[3];
1938 
1939  // calculate new cluster position
1940 
1941  recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1942  recPoint->GetGlobalPosition(gpos);
1943  gpos.GetXYZ(g);
1944 
1945  // create a new cluster
1946 
1947  (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1948  AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1949  j++;
1950  clus->SetType(AliVCluster::kEMCALClusterv1);
1951  clus->SetE(clusterE);
1952  clus->SetPosition(g);
1953  clus->SetNCells(ncellsTrue);
1954  clus->SetCellsAbsId(absIds);
1955  clus->SetCellsAmplitudeFraction(ratios);
1956  clus->SetChi2(-1); //not yet implemented
1957  clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1958  clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1959  clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1960 
1961  if(ncells == ncellsTrue)
1962  {
1963  Float_t elipAxis[2];
1964  recPoint->GetElipsAxis(elipAxis);
1965  clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1966  clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1967  clus->SetDispersion(recPoint->GetDispersion());
1968  }
1969  else if(fSelectCell)
1970  {
1971  // In case some cells rejected, in unfolding case, recalculate
1972  // shower shape parameters and position
1973  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
1974 
1975  AliVCaloCells* cells = 0x0;
1976  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1977  else cells = InputEvent()->GetEMCALCells();
1978 
1979  fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1980  fRecoUtils->RecalculateClusterPID(clus);
1981  fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1982  }
1983 
1984  // MC
1985 
1988  else
1989  {
1990  //
1991  // Normal case, trust what the clusterizer has found
1992  //
1993  Int_t parentMult = 0;
1994  Int_t *parentList = recPoint->GetParents(parentMult);
1995  Float_t *parentListDE = recPoint->GetParentsDE(); // deposited energy
1996 
1997  clus->SetLabel(parentList, parentMult);
1999  clus->SetClusterMCEdepFractionFromEdepArray(parentListDE);
2000 
2001  //
2002  // Set the cell energy deposition fraction map:
2003  //
2004  if( parentMult > 0 && fSetCellMCLabelFromEdepFrac )
2005  {
2006  UInt_t * mcEdepFracPerCell = new UInt_t[ncellsTrue];
2007 
2008  // Get the digit that originated this cell cluster
2009 // AliVCaloCells* cells = 0x0;
2010 // if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
2011 // else cells = InputEvent()->GetEMCALCells();
2012 
2013  for(Int_t icell = 0; icell < ncellsTrue ; icell++)
2014  {
2015  Int_t idigit = fCellLabels[absIds[icell]];
2016 
2017  const AliEMCALDigit * dig = (const AliEMCALDigit*)fDigitsArr->At(idigit);
2018 
2019  // Find the 4 MC labels that contributed to the cluster and their
2020  // deposited energy in the current digit
2021 
2022  mcEdepFracPerCell[icell] = 0; // init
2023 
2024  Int_t nparents = dig->GetNiparent();
2025  if ( nparents > 0 )
2026  {
2027  Int_t digLabel =-1 ;
2028  Float_t edep = 0 ;
2029  Float_t edepTot = 0 ;
2030  Float_t mcEDepFrac[4] = {0,0,0,0};
2031 
2032  // all parents in digit
2033  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2034  {
2035  digLabel = dig->GetIparent (jndex+1);
2036  edep = dig->GetDEParent(jndex+1);
2037  edepTot += edep;
2038 
2039  if ( digLabel == parentList[0] ) mcEDepFrac[0] = edep;
2040  else if ( digLabel == parentList[1] ) mcEDepFrac[1] = edep;
2041  else if ( digLabel == parentList[2] ) mcEDepFrac[2] = edep;
2042  else if ( digLabel == parentList[3] ) mcEDepFrac[3] = edep;
2043  } // all prarents in digit
2044 
2045  // Divide energy deposit by total deposited energy
2046  // Do this only when deposited energy is significant, use 10 MeV although 50 MeV should be expected
2047  if(edepTot > 0.01)
2048  {
2049  mcEdepFracPerCell[icell] = clus->PackMCEdepFraction(mcEDepFrac);
2050  }
2051  } // at least one parent label in digit
2052  } // cell in cluster loop
2053 
2054  clus->SetCellsMCEdepFractionMap(mcEdepFracPerCell);
2055 
2056  delete [] mcEdepFracPerCell;
2057 
2058  } // at least one parent in cluster, do the cell primary packing
2059  }
2060  } // recPoints loop
2061 }
2062 
2063 //_____________________________________________________________________
2066 //_____________________________________________________________________
2068 {
2069  if(label < 0) return ;
2070 
2071  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
2072  if(!evt) return ;
2073 
2074  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
2075  if(!arr) return ;
2076 
2077  if(label < arr->GetEntriesFast())
2078  {
2079  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
2080  if(!particle) return ;
2081 
2082  if(label == particle->Label()) return ; // label already OK
2083  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
2084  }
2085  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
2086 
2087  // loop on the particles list and check if there is one with the same label
2088  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
2089  {
2090  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
2091  if(!particle) continue ;
2092 
2093  if(label == particle->Label())
2094  {
2095  label = ind;
2096  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
2097  return;
2098  }
2099  }
2100 
2101  label = -1;
2102 
2103  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
2104 }
2105 
2106 //________________________________________________
2108 //________________________________________________
2110 {
2111  for(Int_t j = 0; j < fgkNEMCalCells; j++)
2112  {
2113  fOrgClusterCellId[j] = -1;
2114  fCellLabels[j] = -1;
2115  fCellSecondLabels[j] = -1;
2116  fCellTime[j] = 0.;
2117  fCellMatchdEta[j] = -999;
2118  fCellMatchdPhi[j] = -999;
2119 
2120  fTCardCorrCellsEner[j] = 0.;
2121  fTCardCorrCellsNew [j] = kFALSE;
2122  }
2123 }
2124 
2125 //_____________________________________________________________________________________________________
2129 //_____________________________________________________________________________________________________
2131  AliAODCaloCluster * clus)
2132 {
2133  Int_t parentMult = 0;
2134  Int_t *parentList = recPoint->GetParents(parentMult);
2135  clus->SetLabel(parentList, parentMult);
2136 
2137  //Write the second major contributor to each MC cluster.
2138  Int_t iNewLabel ;
2139  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2140  {
2141 
2142  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2143  if(idCell>=0)
2144  {
2145  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
2146  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
2147  {
2148  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
2149  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
2150  }
2151  if (iNewLabel == 1)
2152  {
2153  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
2154  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
2155  {
2156  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
2157  }
2158 
2159  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
2160  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
2161  delete [] newLabelArray;
2162  }
2163  } // positive cell number
2164  }
2165 }
2166 
2167 //___________________________________________________________________________________________________
2171 //___________________________________________________________________________________________________
2173 {
2174  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
2175  clArray.Reset();
2176  Int_t nClu = 0;
2177  Int_t nLabTotOrg = 0;
2178  Float_t emax = -1;
2179  Int_t idMax = -1;
2180 
2181  AliVEvent * event = fEvent;
2182 
2183  // In case of embedding the MC signal is in the added event
2184  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2185  if(aodIH && aodIH->GetEventToMerge()) //Embedding
2186  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
2187 
2188 
2189  // Find the clusters that originally had the cells
2190  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
2191  {
2192  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
2193 
2194  if(idCell>=0)
2195  {
2196  Int_t idCluster = fOrgClusterCellId[idCell];
2197 
2198  Bool_t set = kTRUE;
2199  for(Int_t icl =0; icl < nClu; icl++)
2200  {
2201  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
2202  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
2203  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
2204  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
2205  }
2206  if( set && idCluster >= 0)
2207  {
2208  clArray.SetAt(idCluster,nClu++);
2209  //printf("******** idCluster %d \n",idCluster);
2210  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
2211 
2212  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
2213 
2214  //Search highest E cluster
2215  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2216  //printf("\t E %f\n",clOrg->E());
2217  if(emax < clOrg->E())
2218  {
2219  emax = clOrg->E();
2220  idMax = idCluster;
2221  }
2222  }
2223  }
2224  }// cell loop
2225 
2226  if(nClu==0 || nLabTotOrg == 0)
2227  {
2228  //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());
2229  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
2230  //printf("\n");
2231  }
2232 
2233  // Put the first in the list the cluster with highest energy
2234  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
2235  {
2236  Int_t maxIndex = -1;
2237  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
2238  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2239  {
2240  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
2241  }
2242 
2243  if(firstCluster >=0 && idMax >=0)
2244  {
2245  clArray.SetAt(idMax,0);
2246  clArray.SetAt(firstCluster,maxIndex);
2247  }
2248  }
2249 
2250  // Get the labels list in the original clusters, assign all to the new cluster
2251  TArrayI clMCArray(nLabTotOrg) ;
2252  clMCArray.Reset();
2253 
2254  Int_t nLabTot = 0;
2255  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
2256  {
2257  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
2258  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
2259  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
2260  Int_t nLab = clOrg->GetNLabels();
2261 
2262  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
2263  {
2264  Int_t lab = clOrg->GetLabelAt(iLab) ;
2265  if(lab>=0)
2266  {
2267  Bool_t set = kTRUE;
2268  //printf("\t \t Set Label %d \n", lab);
2269  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
2270  {
2271  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
2272  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
2273  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
2274  }
2275  if( set ) clMCArray.SetAt(lab,nLabTot++);
2276  }
2277  }
2278  }// cluster loop
2279 
2280  // Set the final list of labels
2281 
2282  clus->SetLabel(clMCArray.GetArray(), nLabTot);
2283 
2284 // printf("Final list of labels for new cluster : \n");
2285 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
2286 // {
2287 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
2288 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
2289 // printf(" org %d ",label);
2290 // RemapMCLabelForAODs(label);
2291 // printf(" new %d \n",label);
2292 // }
2293 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
2294 }
2295 
2296 //____________________________________________________________
2297 // Init geometry, create list of output clusters.
2298 //____________________________________________________________
2300 {
2301  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2302 
2303  if(fOutputAODBranchName.Length()==0)
2304  {
2305  fOutputAODBranchName = "newEMCALClustersArray";
2306  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
2307  }
2308 
2310 
2311  if( fFillAODFile )
2312  {
2313  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
2314 
2315  //fOutputAODBranch->SetOwner(kFALSE);
2316 
2317  AddAODBranch("TClonesArray", &fOutputAODBranch);
2318  }
2319 }
2320 
2321 //_______________________________________________________________________
2324 //________________________________________________________________________
2326 {
2327  if ( !fUpdateCell ) return;
2328 
2329  // Update cells only in case re-calibration was done
2330  // or bad map applied or additional T-Card cells added.
2331  if(!fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
2332  !fRecoUtils->IsRecalibrationOn() &&
2333  !fRecoUtils->IsRunDepRecalibrationOn() &&
2334  !fRecoUtils->IsTimeRecalibrationOn() &&
2335  !fRecoUtils->IsL1PhaseInTimeRecalibrationOn() &&
2336  !fTCardCorrEmulation ) return;
2337 
2338  const Int_t ncells = fCaloCells->GetNumberOfCells();
2339  const Int_t ndigis = fDigitsArr->GetEntries();
2340  if ( ncells != ndigis )
2341  {
2342  fCaloCells->DeleteContainer();
2343  fCaloCells->CreateContainer(ndigis);
2344  }
2345 
2346  for (Int_t idigit = 0; idigit < ndigis; ++idigit)
2347  {
2348  AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
2349 
2350  Double_t cellAmplitude = digit->GetAmplitude();
2351  Short_t cellNumber = digit->GetId();
2352  Double_t cellTime = digit->GetTime();
2353 
2354  Bool_t highGain = kFALSE;
2355  if( digit->GetType() == AliEMCALDigit::kHG ) highGain = kTRUE;
2356 
2357  // Only for MC
2358  // Get the label of the primary particle that generated the cell
2359  // Assign the particle that deposited more energy
2360  Int_t nparents = digit->GetNiparent();
2361  Int_t cellMcEDepFrac =-1 ;
2362  Float_t cellMcLabel =-1.;
2363  if ( nparents > 0 )
2364  {
2365  for ( Int_t jndex = 0 ; jndex < nparents ; jndex++ )
2366  {
2367  if(cellMcEDepFrac >= digit->GetDEParent(jndex+1)) continue ;
2368 
2369  cellMcLabel = digit->GetIparent (jndex+1);
2370  cellMcEDepFrac= digit->GetDEParent(jndex+1);
2371  } // all primaries in digit
2372  } // select primary label
2373 
2374  if ( cellMcEDepFrac < 0 ) cellMcEDepFrac = 0.;
2375 
2376  fCaloCells->SetCell(idigit, cellNumber, cellAmplitude, cellTime, cellMcLabel, cellMcEDepFrac, highGain);
2377  }
2378 }
2379 
2380 
2381 //_______________________________________________________
2392 //_______________________________________________________
2394 {
2395  //-------
2396  // Step 1
2397 
2398  // Remove the contents of AOD branch output list set in the previous event
2399  fOutputAODBranch->Clear("C");
2400 
2401  LoadBranches();
2402 
2403  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
2404  // If we need a centrality bin, we select only those events in the corresponding bin.
2405  if( fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
2406  {
2407  Float_t cen = GetEventCentrality();
2408  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
2409  }
2410 
2411  // Intermediate array with new clusters : init the array only once or clear from previous event
2412  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
2413  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
2414 
2415  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
2416 
2417  // Get the event, do some checks and settings
2418  CheckAndGetEvent() ;
2419 
2420  if (!fEvent)
2421  {
2422  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
2423  return ;
2424  }
2425 
2426  // Init pointers, geometry, clusterizer, ocdb, aodb
2427 
2428  if(fAccessOCDB) AccessOCDB();
2429  if(fAccessOADB) AccessOADB(); // only once
2430 
2432 
2433  //-------
2434  // Step 2
2435 
2436  // Make clusters
2438  else ClusterizeCells() ;
2439 
2440  //-------
2441  // Step 3
2442 
2444 
2445  UpdateCells();
2446 }
2447 
2448 
void FillAODCaloCells()
Put calo cells in standard branch.
Float_t fTCardCorrInduceEnerFracWidth[4]
Induced energy loss gauss witdth on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ ...
AliVCaloCells * fCaloCells
! CaloCells container
Bool_t fTCardCorrEmulation
Activate T-Card cells energy correlation.
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
Bool_t fTCardCorrClusEnerConserv
When making correlation, subtract from the reference cell the induced energy on the neighbour cells...
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.
Float_t fTCardCorrCellsEner[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells.
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.
Float_t fTCardCorrInduceEnerFrac[4]
Induced energy loss gauss mean on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ co...
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 fTCardCorrInduceEnerFracP1[4]
Induced energy loss gauss mean on 0-same row, diff col, 1-up/down cells left/right col 2-left/righ co...
Float_t fConstantTimeShift
Apply a 600 ns time shift in case of simulation, shift in ns.
void SetClustersMCLabelFromOriginalClusters(AliAODCaloCluster *clus)
Bool_t fUpdateCell
On/Off the upate of the CaloCells container.
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 recalculation to speed up analysis in PbPb.
Float_t fCellMatchdPhi[fgkNEMCalCells]
Array with cluster-track dEta.
Float_t fTCardCorrInduceEnerProb[22]
Probability to induce energy loss per SM.
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
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
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.
TFile * file
TList with histograms for a given trigger.
Bool_t fAccessOCDB
Need to access info from OCDB (not really)
TObjArray * fClusterArr
! Recpoints array
unsigned short UShort_t
Definition: External.C:28
Bool_t fRecalibrateWithClusterTime
Use fCellTime to store time of cells in cluster.
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
void FillAODHeader()
Put event header information in standard AOD branch.
static const Int_t fgkNEMCalCells
Total number of cells in the calorimeter, 10*48*24 (EMCal) + 4*48*8 (EMCal/DCal 1/3) + 6*32*24 (DCal)...
Int_t fEMCALNcellsCut
At least an EMCAL cluster with fNCellsCut cells over fEnergyCut.
Reclusterize EMCal clusters, put them in a new branch for other following analysis.
Bool_t fTCardCorrCellsNew[fgkNEMCalCells]
Array with induced cell energy in T-Card neighbour cells, that before had no signal.
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.