AliPhysics  a0db429 (a0db429)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
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 
42 // --- EMCAL
43 #include "AliEMCALAfterBurnerUF.h"
44 #include "AliEMCALGeometry.h"
45 #include "AliEMCALClusterizerNxN.h"
46 #include "AliEMCALClusterizerv1.h"
47 #include "AliEMCALClusterizerv2.h"
48 #include "AliEMCALRecPoint.h"
49 #include "AliEMCALDigit.h"
50 
52 
56 
57 //______________________________________________________________________________
58 // Constructor.
59 //______________________________________________________________________________
61 : AliAnalysisTaskSE(name)
62 , fEvent(0)
63 , fGeom(0), fGeomName("")
64 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
65 , fOCDBpath(""), fAccessOCDB(kFALSE)
66 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
67 , fRecParam(0), fClusterizer(0)
68 , fUnfolder(0), fJustUnfold(kFALSE)
69 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
70 , fFillAODFile(kFALSE), fFillAODHeader(0)
71 , fFillAODCaloCells(0), fRun(-1)
72 , fRecoUtils(0), fConfigName("")
73 , fOrgClusterCellId()
74 , fCellLabels(), fCellSecondLabels(), fCellTime()
75 , fCellMatchdEta(), fCellMatchdPhi()
76 , fRecalibrateWithClusterTime(0)
77 , fMaxEvent(0), fDoTrackMatching(kFALSE)
78 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
79 , fRejectBelowThreshold(kFALSE)
80 , fRemoveLEDEvents(kTRUE),fRemoveExoticEvents(kFALSE)
81 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
82 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
83 , fCentralityClass(""), fSelectEMCALEvent(0)
84 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
85 , fSetCellMCLabelFromCluster(0)
86 , fRemapMCLabelForAODs(0)
87 , fInputFromFilter(0)
88 {
89  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = 0;
90 
91  ResetArrays();
92 
93  fCentralityBin[0] = fCentralityBin[1]=-1;
94 }
95 
96 //______________________________________________________________
98 //______________________________________________________________
100 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
101 , fEvent(0)
102 , fGeom(0), fGeomName("")
103 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
104 , fOCDBpath(""), fAccessOCDB(kFALSE)
105 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
106 , fRecParam(0), fClusterizer(0)
107 , fUnfolder(0), fJustUnfold(kFALSE)
108 , fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
109 , fFillAODFile(kFALSE), fFillAODHeader(0)
110 , fFillAODCaloCells(0), fRun(-1)
111 , fRecoUtils(0), fConfigName("")
112 , fOrgClusterCellId()
113 , fCellLabels(), fCellSecondLabels(), fCellTime()
114 , fCellMatchdEta(), fCellMatchdPhi()
115 , fRecalibrateWithClusterTime(0)
116 , fMaxEvent(0), fDoTrackMatching(kFALSE)
117 , fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
118 , fRejectBelowThreshold(kFALSE)
119 , fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
120 , fImportGeometryFromFile(kTRUE), fImportGeometryFilePath("")
121 , fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
122 , fCentralityClass(""), fSelectEMCALEvent(0)
123 , fEMCALEnergyCut(0.), fEMCALNcellsCut (0)
124 , fSetCellMCLabelFromCluster(0)
125 , fRemapMCLabelForAODs(0)
126 , fInputFromFilter(0)
127 {
128  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = 0;
129 
130  ResetArrays();
131 
132  fCentralityBin[0] = fCentralityBin[1]=-1;
133 }
134 
135 //_______________________________________________________________
137 //_______________________________________________________________
139 {
140  if (fDigitsArr)
141  {
142  fDigitsArr->Clear("C");
143  delete fDigitsArr;
144  }
145 
146  if (fClusterArr)
147  {
148  fClusterArr->Delete();
149  delete fClusterArr;
150  }
151 
152  if (fCaloClusterArr)
153  {
154  fCaloClusterArr->Delete();
155  delete fCaloClusterArr;
156  }
157 
158  if(fClusterizer) delete fClusterizer;
159  if(fUnfolder) delete fUnfolder;
160  if(fRecoUtils) delete fRecoUtils;
161 }
162 
163 //_______________________________________________________
167 //_______________________________________________________
169 {
170  if(!fSelectEMCALEvent) return kTRUE; // accept
171 
172  if(fEMCALEnergyCut <= 0) return kTRUE; // accept
173 
174  Int_t nCluster = fEvent -> GetNumberOfCaloClusters();
175  AliVCaloCells * caloCell = fEvent -> GetEMCALCells();
176  Int_t bc = fEvent -> GetBunchCrossNumber();
177 
178  for(Int_t icalo = 0; icalo < nCluster; icalo++)
179  {
180  AliVCluster *clus = (AliVCluster*) (fEvent->GetCaloCluster(icalo));
181 
182  if( ( clus->IsEMCAL() ) && ( clus->GetNCells() > fEMCALNcellsCut ) && ( clus->E() > fEMCALEnergyCut ) &&
183  fRecoUtils->IsGoodCluster(clus,fGeom,caloCell,bc))
184  {
185 
186  AliDebug(1, Form("Accept : E %2.2f > %2.2f, nCells %d > %d",
187  clus->E(), fEMCALEnergyCut, clus->GetNCells(), fEMCALNcellsCut));
188 
189  return kTRUE;
190  }
191 
192  }// loop
193 
194  AliDebug(1,"Reject");
195 
196  return kFALSE;
197 }
198 
199 //_______________________________________________
202 //_______________________________________________
204 {
205  //Set it only once
206  if(fOADBSet) return ;
207 
208  Int_t runnumber = InputEvent()->GetRunNumber() ;
209  TString pass = GetPass();
210 
211  AliInfo(Form("Get AODB parameters from EMCAL in %s for run %d, and <%s>",fOADBFilePath.Data(),runnumber,pass.Data()));
212 
213  Int_t nSM = fGeom->GetNumberOfSuperModules();
214 
215  // Bad map
216  if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
217  {
218  AliOADBContainer *contBC=new AliOADBContainer("");
219  contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
220 
221  TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
222 
223  if(arrayBC)
224  {
225  fRecoUtils->SwitchOnDistToBadChannelRecalculation();
226  AliInfo("Remove EMCAL bad cells");
227 
228  for (Int_t i=0; i < nSM; ++i)
229  {
230  TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
231 
232  if (hbm)
233  delete hbm;
234 
235  hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
236 
237  if (!hbm)
238  {
239  AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
240  continue;
241  }
242 
243  hbm->SetDirectory(0);
244  fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
245 
246  } // loop
247  } else AliInfo("Do NOT remove EMCAL bad channels"); // run array
248 
249  delete contBC;
250  } // Remove bad
251 
252  // Energy Recalibration
253  if(fRecoUtils->IsRecalibrationOn())
254  {
255  AliOADBContainer *contRF=new AliOADBContainer("");
256 
257  contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
258 
259  TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
260 
261  if(recal)
262  {
263  TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
264 
265  if(recalpass)
266  {
267  TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
268 
269  if(recalib)
270  {
271  AliInfo("Recalibrate EMCAL");
272  for (Int_t i=0; i<nSM; ++i)
273  {
274  TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
275 
276  if (h)
277  delete h;
278 
279  h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
280 
281  if (!h)
282  {
283  AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
284  continue;
285  }
286 
287  h->SetDirectory(0);
288 
289  fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
290  } // SM loop
291  } else AliInfo("Do NOT recalibrate EMCAL, no params object array"); // array ok
292  } else AliInfo("Do NOT recalibrate EMCAL, no params for pass"); // array pass ok
293  } else AliInfo("Do NOT recalibrate EMCAL, no params for run"); // run number array ok
294 
295  delete contRF;
296  } // Recalibration on
297 
298  // Energy Recalibration, apply on top of previous calibration factors
299  if(fRecoUtils->IsRunDepRecalibrationOn())
300  {
301  AliOADBContainer *contRFTD=new AliOADBContainer("");
302 
303  contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
304 
305  TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
306 
307  //If it did not exist for this run, get closes one
308  if (!htd)
309  {
310  AliWarning(Form("No TemperatureCorrCalib Objects for run: %d",runnumber));
311  // let's get the closest runnumber instead then..
312  Int_t lower = 0;
313  Int_t ic = 0;
314  Int_t maxEntry = contRFTD->GetNumberOfEntries();
315 
316  while ( (ic < maxEntry) && (contRFTD->UpperLimit(ic) < runnumber) ) {
317  lower = ic;
318  ic++;
319  }
320 
321  Int_t closest = lower;
322  if ( (ic<maxEntry) &&
323  (contRFTD->LowerLimit(ic)-runnumber) < (runnumber - contRFTD->UpperLimit(lower)) ) {
324  closest = ic;
325  }
326 
327  AliWarning(Form("TemperatureCorrCalib Objects found closest id %d from run: %d", closest, contRFTD->LowerLimit(closest)));
328  htd = (TH1S*) contRFTD->GetObjectByIndex(closest);
329  }
330 
331  if(htd)
332  {
333  AliInfo("Recalibrate (Temperature) EMCAL");
334 
335  for (Int_t ism=0; ism<nSM; ++ism)
336  {
337  for (Int_t icol=0; icol<48; ++icol)
338  {
339  for (Int_t irow=0; irow<24; ++irow)
340  {
341  Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
342 
343  Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
344  factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
345  //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
346  // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
347  fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
348  } // columns
349  } // rows
350  } // SM loop
351  } else AliInfo("Do NOT recalibrate EMCAL with T variations, no params TH1");
352 
353  delete contRFTD;
354  } // Run by Run T calibration
355 
356  // Time Recalibration
357  if(fRecoUtils->IsTimeRecalibrationOn())
358  {
359  AliOADBContainer *contTRF=new AliOADBContainer("");
360 
361  contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
362 
363  TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
364 
365  if(trecal)
366  {
367  TString passM = pass;
368  if(pass=="spc_calo") passM = "pass1";
369  TObjArray *trecalpass=(TObjArray*)trecal->FindObject(passM);
370 
371  if(trecalpass)
372  {
373  AliInfo("Time Recalibrate EMCAL");
374  for (Int_t ibc = 0; ibc < 4; ++ibc)
375  {
376  TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
377 
378  if (h)
379  delete h;
380 
381  h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
382 
383  if (!h)
384  {
385  AliError(Form("Could not load hAllTimeAvBC%d",ibc));
386  continue;
387  }
388 
389  h->SetDirectory(0);
390 
391  fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
392  } // bunch crossing loop
393  } else AliInfo("Do NOT recalibrate time EMCAL, no params for pass"); // array pass ok
394  } else AliInfo("Do NOT recalibrate time EMCAL, no params for run"); // run number array ok
395 
396  delete contTRF;
397  } // Time recalibration on
398 
399  // Parameters already set once, so do not it again
400  fOADBSet = kTRUE;
401 }
402 
403 //_________________________________________________
406 //_________________________________________________
408 {
409  fEvent = InputEvent();
410  if (!fEvent)
411  {
412  Warning("AccessOCDB","Event not available!!!");
413  return kFALSE;
414  }
415 
416  if (fEvent->GetRunNumber()==fRun)
417  return kTRUE;
418  fRun = fEvent->GetRunNumber();
419 
420  AliDebug(1,"Begin");
421 
422  AliCDBManager *cdb = AliCDBManager::Instance();
423 
424 
425  if (fOCDBpath.Length())
426  {
427  cdb->SetDefaultStorage(fOCDBpath.Data());
428  AliInfo(Form("Default storage %s",fOCDBpath.Data()));
429  }
430 
431  cdb->SetRun(fEvent->GetRunNumber());
432 
433  //
434  // EMCAL from RAW OCDB
435  if (fOCDBpath.Contains("alien:"))
436  {
437  cdb->SetSpecificStorage("EMCAL/Calib/Data","raw://");
438  cdb->SetSpecificStorage("EMCAL/Calib/Time","raw://");
439  cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","raw://");
440  }
441 
442  TString path = cdb->GetDefaultStorage()->GetBaseFolder();
443 
444  return kTRUE;
445 }
446 
447 //_____________________________________________________
452 //_____________________________________________________
454 {
455  fEvent = 0x0;
456 
457  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
458  Int_t eventN = Entry();
459  if(aodIH) eventN = aodIH->GetReadEntry();
460 
461  if (eventN > fMaxEvent)
462  return ;
463 
464  //printf("Clusterizer --- Event %d-- \n",eventN);
465 
466  // Check if input event are embedded events
467  // If so, take output event
468  if (aodIH && aodIH->GetMergeEvents())
469  {
470  fEvent = AODEvent();
471 
472  if(!aodIH->GetMergeEMCALCells())
473  AliFatal("Events merged but not EMCAL cells, check analysis settings!");
474 
475  AliDebug(1,"Use embedded events");
476 
477  AliDebug(1,Form("\t InputEvent N Clusters %d, N Cells %d",
478  InputEvent()->GetNumberOfCaloClusters(),InputEvent()->GetEMCALCells()->GetNumberOfCells()));
479 
480  AliDebug(1,Form("\t MergedEvent N Clusters %d, N Cells %d",
481  aodIH->GetEventToMerge()->GetNumberOfCaloClusters(), aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells()));
482 
483 // if(DebugLevel() > 1)
484 // {
485 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
486 // {
487 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
488 // if(sigCluster->IsEMCAL()) AliInfo(Form("\t \t Signal cluster: i %d, E %f",icl,sigCluster->E()));
489 // }
490 // }
491 
492  AliDebug(1,Form("\t OutputEvent N Clusters %d, N Cells %d",
493  AODEvent()->GetNumberOfCaloClusters(), AODEvent()->GetEMCALCells()->GetNumberOfCells()));
494  }
495  else if(fInputFromFilter)
496  {
497  //printf("Get Input From Filtered AOD\n");
498  fEvent = AODEvent();
499  }
500  else
501  {
502  fEvent = InputEvent();
505  }
506 
507  if (!fEvent)
508  {
509  AliError("Event not available");
510  return ;
511  }
512 
513  //Process events if there is a high energy cluster
514  if(!AcceptEventEMCAL()) { fEvent = 0x0 ; return ; }
515 
516  //-------------------------------------------------------------------------------------
517  // Reject events if LED was firing, use only for LHC11a data
518  // Reject event if triggered by exotic cell and remove exotic cells if not triggered
519  //-------------------------------------------------------------------------------------
520 
521  if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
522 
523  if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
524 
525  //-------------------------------------------------------------------------------------
526  // Set the cluster array in the event (output or input)
527  //-------------------------------------------------------------------------------------
528 
529  if ( fFillAODFile )
530  {
531  //Magic line to write events to AOD filem put after event rejection
532  AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
533  }
534  else if( !fOutputAODBranchSet )
535  {
536  // Create array and put it in the input event, if output AOD not selected, only once
537  InputEvent()->AddObject(fOutputAODBranch);
538  fOutputAODBranchSet = kTRUE;
539  AliInfo(Form("Add AOD branch <%s> to input event",fOutputAODBranchName.Data()));
540  }
541 }
542 
543 //____________________________________________________
548 //____________________________________________________
550 {
551  Int_t nClusters = fEvent->GetNumberOfCaloClusters();
552  Int_t nClustersOrg = 0;
553 
554  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
555  if(aodIH && aodIH->GetEventToMerge()) //Embedding
556  nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
557 
558  ResetArrays();
559 
560  // Loop on original clusters, get MC labels.
561  for (Int_t i = 0; i < nClusters; i++)
562  {
563  AliVCluster *clus = 0;
564  if(aodIH && aodIH->GetEventToMerge()) //Embedding
565  clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
566  else
567  clus = fEvent->GetCaloCluster(i);
568 
569  if(!clus) return;
570 
571  if(clus->IsEMCAL())
572  {
573  Int_t label = clus->GetLabel();
574  Int_t label2 = -1 ;
575  //printf("Org cluster E %f, Time %e, Index = %d, ID %d, MC label %d\n", clus->E(), clus->GetTOF(),i, clus->GetID(),label );
576  //printf("Original list of labels from old cluster : \n");
577  //for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
578 
579  if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
580  UShort_t * index = clus->GetCellsAbsId() ;
581  for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
582  {
583  //printf("\t cell %d, MC label %d\n",index[icell],fEvent->GetEMCALCells()->GetCellMCLabel(index[icell]));
584  fOrgClusterCellId[index[icell]] = i;
585  fCellLabels[index[icell]] = label;
586  fCellSecondLabels[index[icell]] = label2;
587  fCellTime[index[icell]] = clus->GetTOF();
588  fCellMatchdEta[index[icell]] = clus->GetTrackDz();
589  fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
590  }
591  nClustersOrg++;
592  }
593  // printf("\n");
594  }
595 
596  // Transform CaloCells into Digits
597 
598  Int_t idigit = 0;
599  Int_t id = -1;
600  Float_t amp = -1;
601  Double_t time = -1;
602 
603  AliVCaloCells * cells = fEvent->GetEMCALCells();
604 
605  Int_t bc = InputEvent()->GetBunchCrossNumber();
606 
607  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
608  {
609  // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
610  id = cells->GetCellNumber(icell);
611  Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
612 
613  // Do not include cells with too low energy, nor exotic cell
614  if( amp < fRecParam->GetMinECut() ||
615  time > fRecParam->GetTimeMax() ||
616  time < fRecParam->GetTimeMin() ) accept = kFALSE;
617 
618  // In case of old AOD analysis cell time is -1 s, approximate replacing by time of the cluster the digit belongs.
620  {
621  time = fCellTime[id];
622  //printf("cell %d time org %f - ",id, time*1.e9);
623  fRecoUtils->RecalibrateCellTime(id,bc,time);
624  //printf("recal %f\n",time*1.e9);
625  }
626 
627  //Exotic?
628  if (accept && fRecoUtils->IsExoticCell(id,cells,bc))
629  accept = kFALSE;
630 
631  if( !accept )
632  {
633  AliDebug(2,Form("Remove channel absId %d, index %d of %d, amp %f, time %f",
634  id,icell, cells->GetNumberOfCells(), amp, time*1.e9));
635  continue;
636  }
637 
638  Int_t mcLabel = cells->GetMCLabel(icell);
639  //printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - cell %d, mc label %d\n",id,mcLabel);
640 
641  //if(fCellLabels[id]!=mcLabel)printf("mcLabel %d - %d\n",mcLabel,fCellLabels[id]);
642  if ( fSetCellMCLabelFromCluster == 1 ) mcLabel = fCellLabels[id]; // Older aliroot MC productions
644  else mcLabel = -1; // found later
645 
646  //printf("\t new label %d\n",mcLabel);
647 
648  // Create the digit, put a fake primary deposited energy to trick the clusterizer
649  // when checking the most likely primary
650 
651  Float_t efrac = cells->GetEFraction(icell);
652 
653  //When checking the MC of digits, give weight to cells with embedded signal
654  if (mcLabel > 0 && efrac < 1.e-6) efrac = 1;
655 
656  //printf("******* Cell %d, id %d, e %f, fraction %f, MC label %d, used MC label %d\n",icell,id,amp,cells->GetEFraction(icell),cells->GetMCLabel(icell),mcLabel);
657 
658  new((*fDigitsArr)[idigit]) AliEMCALDigit( mcLabel, mcLabel, id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, amp*efrac);
659  // Last parameter should be MC deposited energy, since it is not available, add just the cell amplitude so that
660  // we give more weight to the MC label of the cell with highest energy in the cluster
661 
662  idigit++;
663  }
664 
665  fDigitsArr->Sort();
666 
667  //-------------------------------------------------------------------------------------
668  //Do the clusterization
669  //-------------------------------------------------------------------------------------
670 
671  fClusterizer->Digits2Clusters("");
672 
673  //-------------------------------------------------------------------------------------
674  //Transform the recpoints into AliVClusters
675  //-------------------------------------------------------------------------------------
676 
678 
679  if(!fCaloClusterArr)
680  {
681  AliWarning(Form("No array with CaloClusters, input RecPoints entries %d",fClusterArr->GetEntriesFast()));
682  return;
683  }
684 
685  AliDebug(1,Form("N clusters: before recluster %d, after recluster %d, recpoints %d",
686  nClustersOrg, fCaloClusterArr->GetEntriesFast(),fClusterArr->GetEntriesFast()));
687 
688 // if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
689 // {
690 // AliInfo("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
691 // fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
692 // fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
693 // }
694 }
695 
696 //_____________________________________________________
698 //_____________________________________________________
700 {
701  AliVCaloCells *cells = fEvent->GetEMCALCells();
702  Double_t cellAmplitude = 0;
703  Double_t cellTime = 0;
704  Short_t cellNumber = 0;
705  Int_t cellMCLabel = 0;
706  Double_t cellEFrac = 0;
707  Int_t nClustersOrg = 0;
708 
709  // Fill the array with the EMCAL clusters, copy them
710  for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
711  {
712  AliVCluster *clus = fEvent->GetCaloCluster(i);
713  if(clus->IsEMCAL())
714  {
715  //recalibrate/remove bad channels/etc if requested
716  if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
717  {
718  continue;
719  }
720 
721  if(fRecoUtils->IsRecalibrationOn())
722  {
723  //Calibrate cluster
724  fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
725 
726  //CalibrateCells
727  for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
728  {
729  if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
730  break;
731 
732  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
733  fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
734  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
735 
736  //Do not include bad channels found in analysis?
737  if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
738  fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
739  continue;
740 
741  cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
742 
743  }// cells loop
744  }// recalibrate
745 
746  //Cast to ESD or AOD, needed to create the cluster array
747  AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
748  AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
749 
750  if (esdCluster)
751  {
752  fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
753  }//ESD
754  else if(aodCluster)
755  {
756  fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
757  }//AOD
758  else
759  AliWarning("Wrong CaloCluster type?");
760 
761  nClustersOrg++;
762  }
763  }
764 
765  //Do the unfolding
766  fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
767 
768  //CLEAN-UP
769  fUnfolder->Clear();
770 }
771 
772 //_____________________________________________________
774 //_____________________________________________________
776 {
777  AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
778  Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
779 
780  AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
781  aodEMcells.CreateContainer(nEMcell);
782  aodEMcells.SetType(AliVCaloCells::kEMCALCell);
783  Double_t calibFactor = 1.;
784  for (Int_t iCell = 0; iCell < nEMcell; iCell++)
785  {
786  Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
787  fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
788  fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
789 
790  if(fRecoUtils->IsRecalibrationOn())
791  {
792  calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
793  }
794 
795  if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
796  { //Channel is not declared as bad
797  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
798  eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
799  }
800  else
801  {
802  aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
803  }
804  }
805 
806  aodEMcells.Sort();
807 }
808 
809 //__________________________________________________
811 //__________________________________________________
813 {
814  AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
815  AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
816 
817  Double_t pos[3] ;
818  Double_t covVtx[6];
819  for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
820 
821  AliAODHeader* header = dynamic_cast<AliAODHeader*>(AODEvent()->GetHeader());
822  if(!header) AliFatal("Not a standard AOD");
823  header->SetRunNumber(fEvent->GetRunNumber());
824 
825  if(esdevent)
826  {
827  TTree* tree = fInputHandler->GetTree();
828  if (tree)
829  {
830  TFile* file = tree->GetCurrentFile();
831  if (file) header->SetESDFileName(file->GetName());
832  }
833  }
834  else if (aodevent) {
835  AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aodevent->GetHeader());
836  if(!aodheader) AliFatal("Not a standard AOD");
837  header->SetESDFileName(aodheader->GetESDFileName());
838  }
839 
840  header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
841  header->SetOrbitNumber(fEvent->GetOrbitNumber());
842  header->SetPeriodNumber(fEvent->GetPeriodNumber());
843  header->SetEventType(fEvent->GetEventType());
844 
845  // Centrality
846  if(fEvent->GetCentrality())
847  header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
848  else
849  header->SetCentrality(0);
850 
851  // Trigger
852  header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
853 
854  header->SetFiredTriggerClasses(fEvent->GetFiredTriggerClasses());
855 
856  header->SetTriggerMask(fEvent->GetTriggerMask());
857  header->SetTriggerCluster(fEvent->GetTriggerCluster());
858 
859  if (esdevent)
860  {
861  header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
862  header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
863  header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
864  }
865  else if (aodevent)
866  {
867  header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
868  header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
869  header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
870  }
871 
872  header->SetMagneticField(fEvent->GetMagneticField());
873  //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
874 
875  header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
876  header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
877  header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
878  header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
879  header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
880 
881  Float_t diamxy[2]={(Float_t)fEvent->GetDiamondX(),(Float_t)fEvent->GetDiamondY()};
882  Float_t diamcov[3];
883  fEvent->GetDiamondCovXY(diamcov);
884  header->SetDiamond(diamxy,diamcov);
885  if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
886  else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
887 
888  //
889  Int_t nVertices = 1 ;/* = prim. vtx*/;
890  Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
891 
892  AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
893 
894  // Access to the AOD container of vertices
895  TClonesArray &vertices = *(AODEvent()->GetVertices());
896  Int_t jVertices=0;
897 
898  // Add primary vertex. The primary tracks will be defined
899  // after the loops on the composite objects (V0, cascades, kinks)
900  fEvent->GetPrimaryVertex()->GetXYZ(pos);
901  Float_t chi = 0;
902  if (esdevent)
903  {
904  esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
905  chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
906  }
907  else if (aodevent)
908  {
909  aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
910  chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
911  }
912 
913  AliAODVertex * primary = new(vertices[jVertices++])
914  AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
915  primary->SetName(fEvent->GetPrimaryVertex()->GetName());
916  primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
917 }
918 
919 //___________________________________________________________
923 //___________________________________________________________
925 {
926  // First recalculate track-matching for the new clusters
927  if(fDoTrackMatching)
928  {
929  fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
930  }
931 
932  // Put the new clusters in the AOD list
933 
934  Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
935 
936  for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
937  {
938  AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
939 
940  newCluster->SetID(i);
941 
942  // Correct cluster energy non linearity
943  newCluster->SetE(fRecoUtils->CorrectClusterEnergyLinearity(newCluster));
944 
945  //Add matched track
946  if(fDoTrackMatching)
947  {
948  Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
949  if(trackIndex >= 0)
950  {
951  newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
952  AliDebug(2,Form("Matched Track index %d to new cluster %d",trackIndex,i));
953  }
954 
955  Float_t dR = 999., dZ = 999.;
956  fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dZ,dR);
957  newCluster->SetTrackDistance(dR,dZ);
958  }
959  else
960  {// Assign previously assigned matched track in reco, very very rough
961  Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
962  newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
963  }
964 
965  //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
966  //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
967  //printf("\n");
968 
969  // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
970  fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
971 
972  new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
973 
974  AliDebug(2,Form("New cluster %d of %d, energy %f, mc label %d",
975  newCluster->GetID(), kNumberOfCaloClusters, newCluster->E(), newCluster->GetLabel()));
976 
977  } // cluster loop
978 
979  fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
980 }
981 
982 //_______________________________________________
984 //_______________________________________________
986 {
987  if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
988  {
989  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null");
990  return TString("");
991  }
992 
993  if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
994  {
995  AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null");
996  return TString("");
997  }
998 
999  TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
1000  if (pass.Contains("ass1")) return TString("pass1");
1001  else if (pass.Contains("ass2")) return TString("pass2");
1002  else if (pass.Contains("ass3")) return TString("pass3");
1003  else if (pass.Contains("ass4")) return TString("pass4");
1004  else if (pass.Contains("ass5")) return TString("pass5");
1005  else if (pass.Contains("LHC11c") && pass.Contains("spc_calo") ) return TString("spc_calo");
1006  else if (pass.Contains("calo") || pass.Contains("high_lumi"))
1007  {
1008  AliInfo("Path contains <calo> or <high-lumi>, set as <pass1>");
1009  return TString("pass1");
1010  }
1011 
1012  // No condition fullfilled, give a default value
1013  AliInfo("Pass number string not found");
1014 
1015  return TString("");
1016 }
1017 
1018 //_________________________________________
1021 //_________________________________________
1023 {
1024  if(fDebug >=0) (AliAnalysisManager::GetAnalysisManager())->AddClassDebug(this->ClassName(),fDebug);
1025 
1026  fOADBSet = kFALSE;
1027  if(fOADBFilePath == "") fOADBFilePath = "$ALICE_PHYSICS/OADB/EMCAL" ;
1028 
1029  fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
1030 
1031  if(!fRecParam) fRecParam = new AliEMCALRecParam;
1032  if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
1033 
1034  if(fMaxEvent <= 0) fMaxEvent = 1000000000;
1035  if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
1036  if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
1037  fRejectBelowThreshold = kFALSE;
1038 
1039  // Centrality
1040  if(fCentralityClass == "") fCentralityClass = "V0M";
1041 
1042  if (fOCDBpath == "") fOCDBpath = "raw://" ;
1043  if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
1044 
1045  if(gROOT->LoadMacro(fConfigName) >=0)
1046  {
1047  AliInfo(Form("Configure analysis with %s",fConfigName.Data()));
1048  AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
1049  fGeomName = clus->fGeomName;
1051  fOCDBpath = clus->fOCDBpath;
1052  fAccessOCDB = clus->fAccessOCDB;
1053  fRecParam = clus->fRecParam;
1054  fJustUnfold = clus->fJustUnfold;
1055  fFillAODFile = clus->fFillAODFile;
1056  fRecoUtils = clus->fRecoUtils;
1057  fConfigName = clus->fConfigName;
1058  fMaxEvent = clus->fMaxEvent;
1061  for(Int_t i = 0; i < 22; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
1063  fCentralityBin[0] = clus->fCentralityBin[0];
1064  fCentralityBin[1] = clus->fCentralityBin[1];
1065  }
1066 }
1067 
1068 //_______________________________________________________
1071 //_______________________________________________________
1073 {
1074  if (fJustUnfold)
1075  {
1076  // init the unfolding afterburner
1077  delete fUnfolder;
1078  fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
1079  return;
1080  }
1081 
1082  //First init the clusterizer
1083  delete fClusterizer;
1084  if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1085  fClusterizer = new AliEMCALClusterizerv1 (fGeom);
1086  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1087  fClusterizer = new AliEMCALClusterizerv2(fGeom);
1088  else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1089  {
1090  fClusterizer = new AliEMCALClusterizerNxN(fGeom);
1091  fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1092  fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
1093  }
1094  else
1095  {
1096  AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1097  }
1098 
1099  //Now set the parameters
1100  fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1101  fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1102  fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1103  fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1104  fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1105  fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1106  fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1107  fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1108  fClusterizer->SetTimeCalibration ( fRecParam->IsTimeCalibrationOn() );
1109  fClusterizer->SetInputCalibrated ( kTRUE );
1110  fClusterizer->SetJustClusters ( kTRUE );
1111 
1112  // Initialize the cluster rec points and digits arrays and get them.
1113  fClusterizer->SetOutput(0);
1114  fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1115  fDigitsArr = fClusterizer->GetDigits();
1116 
1117  //In case of unfolding after clusterization is requested, set the corresponding parameters
1118  if(fRecParam->GetUnfold())
1119  {
1120  Int_t i=0;
1121  for (i = 0; i < 8; i++)
1122  {
1123  fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1124  }//end of loop over parameters
1125 
1126  for (i = 0; i < 3; i++)
1127  {
1128  fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1129  fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1130  }//end of loop over parameters
1131  fClusterizer->SetRejectBelowThreshold(fRejectBelowThreshold);//here we set option of unfolding: split or reject energy
1132  fClusterizer->InitClusterUnfolding();
1133 
1134  }// to unfold
1135 }
1136 
1137 //_________________________________________________
1141 //_________________________________________________
1143 {
1144  if(fGeomMatrixSet) return;
1145 
1146  Int_t runnumber = InputEvent()->GetRunNumber() ;
1147  if (!fGeom)
1148  {
1149  if(fGeomName=="")
1150  {
1151  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(runnumber);
1152  AliInfo(Form("Get EMCAL geometry name <%s> for run %d",fGeom->GetName(),runnumber));
1153  }
1154  else
1155  {
1156  fGeom = AliEMCALGeometry::GetInstance(fGeomName);
1157  AliInfo(Form("Set EMCAL geometry name to <%s>",fGeomName.Data()));
1158  }
1159 
1160  // Init geometry, I do not like much to do it like this ...
1161  if(fImportGeometryFromFile && !gGeoManager)
1162  {
1163  if(fImportGeometryFilePath=="") // If not specified, set location depending on run number
1164  {
1165  // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
1166  if (runnumber < 140000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2010.root";
1167  else if (runnumber < 171000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2011.root";
1168  else if (runnumber < 198000) fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2012.root"; // 2012-2013
1169  else fImportGeometryFilePath = "$ALICE_PHYSICS/OADB/EMCAL/geometry_2015.root"; // >=2015
1170  }
1171 
1172  AliInfo(Form("Import %s",fImportGeometryFilePath.Data()));
1173 
1174  TGeoManager::Import(fImportGeometryFilePath) ;
1175  }
1176 
1177  AliDebug(1,Form("Init for run=%d",runnumber));
1178  if (!gGeoManager) AliDebug(1,"Careful!, gGeoManager not loaded, load misalign matrices");
1179  } // geometry pointer did not exist before
1180 
1181  if(fLoadGeomMatrices)
1182  {
1183  AliInfo("Load user defined EMCAL geometry matrices");
1184 
1185  // OADB if available
1186  AliOADBContainer emcGeoMat("AliEMCALgeo");
1187  emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1188  TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1189 
1190  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1191  {
1192 
1193  if (!fGeomMatrix[mod]) // Get it from OADB
1194  {
1195  AliDebug(2,Form("EMCAL matrices SM %d, %p",mod,((TGeoHMatrix*) matEMCAL->At(mod))));
1196  //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1197 
1198  fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1199  }
1200 
1201  if(fGeomMatrix[mod])
1202  {
1203  if(DebugLevel() > 1)
1204  fGeomMatrix[mod]->Print();
1205 
1206  fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1207  }
1208  else if(gGeoManager)
1209  {
1210  AliWarning(Form("Set matrix for SM %d from gGeoManager",mod));
1211  fGeom->SetMisalMatrix(fGeom->GetMatrixForSuperModuleFromGeoManager(mod),mod) ;
1212  }
1213  else
1214  {
1215  AliError(Form("Alignment atrix for SM %d is not available",mod));
1216  }
1217 
1218  fGeomMatrixSet=kTRUE;
1219 
1220  }//SM loop
1221  }//Load matrices
1222  else if(!gGeoManager)
1223  {
1224  AliInfo("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1225  //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1226  if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1227  {
1228  AliWarning("Use ideal geometry, values geometry matrix not kept in AODs");
1229  }//AOD
1230  else
1231  {
1232  AliDebug(1,"Load Misaligned matrices");
1233 
1234  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1235 
1236  if(!esd)
1237  {
1238  AliError("This event does not contain ESDs?");
1239  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
1240  return;
1241  }
1242 
1243  for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1244  {
1245  if(DebugLevel() > 1)
1246  esd->GetEMCALMatrix(mod)->Print();
1247 
1248  if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1249 
1250  }
1251 
1252  fGeomMatrixSet=kTRUE;
1253 
1254  } // ESD
1255  } // Load matrices from Data
1256 }
1257 
1258 //____________________________________________________
1263 //____________________________________________________
1265 {
1266  if(!fRemoveExoticEvents) return kFALSE;
1267 
1268  // Loop on cells
1269 
1270  AliVCaloCells * cells = fEvent->GetEMCALCells();
1271  Float_t totCellE = 0;
1272  Int_t bc = InputEvent()->GetBunchCrossNumber();
1273 
1274  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1275  {
1276  Float_t ecell = 0 ;
1277  Double_t tcell = 0 ;
1278 
1279  Int_t absID = cells->GetCellNumber(icell);
1280  Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1281  if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
1282  }
1283 
1284  // TString triggerclasses = event->GetFiredTriggerClasses();
1285  // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
1286  // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1287  // return;
1288  //
1289 
1290  //printf("TotE cell %f\n",totCellE);
1291  if(totCellE < 1) return kTRUE;
1292 
1293  return kFALSE;
1294 }
1295 
1296 //________________________________________________________________
1298 //________________________________________________________________
1300 {
1301  if(!fRemoveLEDEvents) return kFALSE;
1302 
1303  //check events of LHC11a period
1304  if(run < 146858 || run > 146860) return kFALSE ;
1305 
1306  // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1307  Int_t ncellsSM3 = 0;
1308  AliVCaloCells * cells = fEvent->GetEMCALCells();
1309  for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1310  {
1311  if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1312  }
1313 
1314  TString triggerclasses = fEvent->GetFiredTriggerClasses();
1315 
1316  Int_t ncellcut = 21;
1317  if(triggerclasses.Contains("EMC")) ncellcut = 35;
1318 
1319  if( ncellsSM3 >= ncellcut)
1320  {
1321  AliInfo(Form("Reject event %d with ncells in SM3 %d",(Int_t)Entry(),ncellsSM3));
1322  if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
1323  return kTRUE;
1324  }
1325 
1326  return kFALSE;
1327 }
1328 
1329 //_______________________________________________________
1333 //_______________________________________________________
1335 {
1336  Int_t j = 0;
1337  for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
1338  {
1339  AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
1340 
1341  const Int_t ncells = recPoint->GetMultiplicity();
1342  Int_t ncellsTrue = 0;
1343 
1344  if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1345 
1346  // cells and their amplitude fractions
1347  UShort_t absIds[ncells];
1348  Double32_t ratios[ncells];
1349 
1350  //For later check embedding
1351  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1352 
1353  Float_t clusterE = 0;
1354  for (Int_t c = 0; c < ncells; c++)
1355  {
1356  AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
1357 
1358  absIds[ncellsTrue] = digit->GetId();
1359  ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
1360 
1361  // In case of unfolding, remove digits with unfolded energy too low
1362  if(fSelectCell)
1363  {
1364  if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1365  {
1366 
1367  AliDebug(2,Form("Too small energy in cell of cluster: cluster cell %f, digit %f",
1368  recPoint->GetEnergiesList()[c],digit->GetAmplitude()));
1369  continue;
1370  } // if cuts
1371  }// Select cells
1372 
1373  //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1374  clusterE +=recPoint->GetEnergiesList()[c];
1375 
1376  // In case of embedding, fill ratio with amount of signal,
1377  if (aodIH && aodIH->GetMergeEvents())
1378  {
1379  //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1380  AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1381  AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1382 
1383  Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1384  //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1385  Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1386  //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1387 
1388  if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1389  //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1390 
1391  }//Embedding
1392 
1393  ncellsTrue++;
1394 
1395  }// cluster cell loop
1396 
1397  if (ncellsTrue < 1)
1398  {
1399  AliDebug(2,Form("Skipping cluster with no cells avobe threshold E = %f, ncells %d",
1400  recPoint->GetEnergy(), ncells));
1401  continue;
1402  }
1403 
1404  //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1405 
1406  if(clusterE < fRecParam->GetClusteringThreshold())
1407  {
1408  AliDebug(2,Form("Remove cluster with energy below seed threshold %f",clusterE));
1409  continue;
1410  }
1411 
1412  TVector3 gpos;
1413  Float_t g[3];
1414 
1415  // calculate new cluster position
1416 
1417  recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
1418  recPoint->GetGlobalPosition(gpos);
1419  gpos.GetXYZ(g);
1420 
1421  // create a new cluster
1422 
1423  (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1424  AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
1425  j++;
1426  clus->SetType(AliVCluster::kEMCALClusterv1);
1427  clus->SetE(clusterE);
1428  clus->SetPosition(g);
1429  clus->SetNCells(ncellsTrue);
1430  clus->SetCellsAbsId(absIds);
1431  clus->SetCellsAmplitudeFraction(ratios);
1432  clus->SetChi2(-1); //not yet implemented
1433  clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1434  clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1435  clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1436 
1437  if(ncells == ncellsTrue)
1438  {
1439  Float_t elipAxis[2];
1440  recPoint->GetElipsAxis(elipAxis);
1441  clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1442  clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1443  clus->SetDispersion(recPoint->GetDispersion());
1444  }
1445  else if(fSelectCell)
1446  {
1447  // In case some cells rejected, in unfolding case, recalculate
1448  // shower shape parameters and position
1449  AliDebug(2,Form("Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape",ncells,ncellsTrue));
1450 
1451  AliVCaloCells* cells = 0x0;
1452  if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1453  else cells = InputEvent()->GetEMCALCells();
1454 
1455  fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1456  fRecoUtils->RecalculateClusterPID(clus);
1457  fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1458 
1459  }
1460 
1461  // MC
1462 
1465  else
1466  {
1467  // Normal case, trust what the clusterizer has found
1468  Int_t parentMult = 0;
1469  Int_t *parentList = recPoint->GetParents(parentMult);
1470  clus->SetLabel(parentList, parentMult);
1471 // printf("Label list : ");
1472 // for(Int_t ilabel = 0; ilabel < parentMult; ilabel++ ) printf(" %d ",parentList[ilabel]);
1473 // printf("\n");
1474  }
1475 
1476  } // recPoints loop
1477 }
1478 
1479 //_____________________________________________________________________
1482 //_____________________________________________________________________
1484 {
1485  if(label < 0) return ;
1486 
1487  AliAODEvent * evt = dynamic_cast<AliAODEvent*> (fEvent) ;
1488  if(!evt) return ;
1489 
1490  TClonesArray * arr = dynamic_cast<TClonesArray*>(evt->FindListObject("mcparticles")) ;
1491  if(!arr) return ;
1492 
1493  if(label < arr->GetEntriesFast())
1494  {
1495  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(label));
1496  if(!particle) return ;
1497 
1498  if(label == particle->Label()) return ; // label already OK
1499  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d - AOD stack %d \n",label, particle->Label());
1500  }
1501  //else printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label %d > AOD labels %d \n",label, arr->GetEntriesFast());
1502 
1503  // loop on the particles list and check if there is one with the same label
1504  for(Int_t ind = 0; ind < arr->GetEntriesFast(); ind++ )
1505  {
1506  AliAODMCParticle * particle = dynamic_cast<AliAODMCParticle *>(arr->At(ind));
1507  if(!particle) continue ;
1508 
1509  if(label == particle->Label())
1510  {
1511  label = ind;
1512  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - New Label Index %d \n",label);
1513  return;
1514  }
1515  }
1516 
1517  label = -1;
1518 
1519  //printf("AliAnalysisTaskEMCALClusterize::RemapMCLabelForAODs() - Label not found set to -1 \n");
1520 }
1521 
1522 //________________________________________________
1524 //________________________________________________
1526 {
1527  for(Int_t j = 0; j < fgkNEMCalCells; j++)
1528  {
1529  fOrgClusterCellId[j] = -1;
1530  fCellLabels[j] = -1;
1531  fCellSecondLabels[j] = -1;
1532  fCellTime[j] = 0.;
1533  fCellMatchdEta[j] = -999;
1534  fCellMatchdPhi[j] = -999;
1535  }
1536 }
1537 
1538 //_____________________________________________________________________________________________________
1542 //_____________________________________________________________________________________________________
1544  AliAODCaloCluster * clus)
1545 {
1546  Int_t parentMult = 0;
1547  Int_t *parentList = recPoint->GetParents(parentMult);
1548  clus->SetLabel(parentList, parentMult);
1549 
1550  //Write the second major contributor to each MC cluster.
1551  Int_t iNewLabel ;
1552  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1553  {
1554 
1555  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1556  if(idCell>=0)
1557  {
1558  iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
1559  for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1560  {
1561  if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
1562  if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
1563  }
1564  if (iNewLabel == 1)
1565  {
1566  Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1567  for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1568  {
1569  newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1570  }
1571 
1572  newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1573  clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1574  delete [] newLabelArray;
1575  }
1576  } // positive cell number
1577  }
1578 }
1579 
1580 //___________________________________________________________________________________________________
1584 //___________________________________________________________________________________________________
1586 {
1587  TArrayI clArray(300) ; //Weird if more than a few clusters are in the origin ...
1588  clArray.Reset();
1589  Int_t nClu = 0;
1590  Int_t nLabTotOrg = 0;
1591  Float_t emax = -1;
1592  Int_t idMax = -1;
1593 
1594  AliVEvent * event = fEvent;
1595 
1596  // In case of embedding the MC signal is in the added event
1597  AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1598  if(aodIH && aodIH->GetEventToMerge()) //Embedding
1599  event = aodIH->GetEventToMerge(); //Get clusters directly from embedded signal
1600 
1601 
1602  // Find the clusters that originally had the cells
1603  for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1604  {
1605  Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1606 
1607  if(idCell>=0)
1608  {
1609  Int_t idCluster = fOrgClusterCellId[idCell];
1610 
1611  Bool_t set = kTRUE;
1612  for(Int_t icl =0; icl < nClu; icl++)
1613  {
1614  if(((Int_t)clArray.GetAt(icl))==-1 ) continue;
1615  if( idCluster == ((Int_t)clArray.GetAt(icl)) ) set = kFALSE;
1616  // printf("\t \t icell %d Cluster in array %d, IdCluster %d, in array %d, set %d\n",
1617  // iLoopCell, icl, idCluster,((Int_t)clArray.GetAt(icl)),set);
1618  }
1619  if( set && idCluster >= 0)
1620  {
1621  clArray.SetAt(idCluster,nClu++);
1622  //printf("******** idCluster %d \n",idCluster);
1623  nLabTotOrg+=(event->GetCaloCluster(idCluster))->GetNLabels();
1624 
1625  //printf("Cluster in array %d, IdCluster %d\n",nClu-1, idCluster);
1626 
1627  //Search highest E cluster
1628  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1629  //printf("\t E %f\n",clOrg->E());
1630  if(emax < clOrg->E())
1631  {
1632  emax = clOrg->E();
1633  idMax = idCluster;
1634  }
1635  }
1636  }
1637  }// cell loop
1638 
1639  if(nClu==0 || nLabTotOrg == 0)
1640  {
1641  //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());
1642  //for(Int_t icell = 0; icell < clus->GetNCells(); icell++) printf("\t cell %d",clus->GetCellsAbsId()[icell]);
1643  //printf("\n");
1644  }
1645 
1646  // Put the first in the list the cluster with highest energy
1647  if(idMax != ((Int_t)clArray.GetAt(0))) // Max not at first position
1648  {
1649  Int_t maxIndex = -1;
1650  Int_t firstCluster = ((Int_t)clArray.GetAt(0));
1651  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1652  {
1653  if(idMax == ((Int_t)clArray.GetAt(iLoopCluster))) maxIndex = iLoopCluster;
1654  }
1655 
1656  if(firstCluster >=0 && idMax >=0)
1657  {
1658  clArray.SetAt(idMax,0);
1659  clArray.SetAt(firstCluster,maxIndex);
1660  }
1661  }
1662 
1663  // Get the labels list in the original clusters, assign all to the new cluster
1664  TArrayI clMCArray(nLabTotOrg) ;
1665  clMCArray.Reset();
1666 
1667  Int_t nLabTot = 0;
1668  for ( Int_t iLoopCluster = 0 ; iLoopCluster < nClu ; iLoopCluster++ )
1669  {
1670  Int_t idCluster = (Int_t) clArray.GetAt(iLoopCluster);
1671  //printf("New Cluster in Array %d, idCluster %d \n",iLoopCluster,idCluster);
1672  AliVCluster * clOrg = event->GetCaloCluster(idCluster);
1673  Int_t nLab = clOrg->GetNLabels();
1674 
1675  for ( Int_t iLab = 0 ; iLab < nLab ; iLab++ )
1676  {
1677  Int_t lab = clOrg->GetLabelAt(iLab) ;
1678  if(lab>=0)
1679  {
1680  Bool_t set = kTRUE;
1681  //printf("\t \t Set Label %d \n", lab);
1682  for(Int_t iLabTot =0; iLabTot < nLabTot; iLabTot++)
1683  {
1684  if( lab == ((Int_t)clMCArray.GetAt(iLabTot)) ) set = kFALSE;
1685  //printf("iLoopCluster %d, Label ID in Org Cluster %d,label %d Label ID in array %d, label in array %d, set %d\n",
1686  // iLoopCluster, iLab, lab, iLabTot, ((Int_t)clMCArray.GetAt(iLabTot)),set);
1687  }
1688  if( set ) clMCArray.SetAt(lab,nLabTot++);
1689  }
1690  }
1691  }// cluster loop
1692 
1693  // Set the final list of labels
1694 
1695  clus->SetLabel(clMCArray.GetArray(), nLabTot);
1696 
1697 // printf("Final list of labels for new cluster : \n");
1698 // for(Int_t ice = 0; ice < clus->GetNCells() ; ice++)
1699 // {
1700 // printf("\t Cell %d ",clus->GetCellsAbsId()[ice]);
1701 // Int_t label = InputEvent()->GetEMCALCells()->GetCellMCLabel(clus->GetCellsAbsId()[ice]);
1702 // printf(" org %d ",label);
1703 // RemapMCLabelForAODs(label);
1704 // printf(" new %d \n",label);
1705 // }
1706 // for(Int_t imc = 0; imc < clus->GetNLabels(); imc++) printf("\t Label %d\n",clus->GetLabelAt(imc));
1707 }
1708 
1709 //____________________________________________________________
1710 // Init geometry, create list of output clusters.
1711 //____________________________________________________________
1713 {
1714  fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1715 
1716  if(fOutputAODBranchName.Length()==0)
1717  {
1718  fOutputAODBranchName = "newEMCALClustersArray";
1719  AliInfo("Cluster branch name not set, set it to newEMCALClustersArray");
1720  }
1721 
1723 
1724  if( fFillAODFile )
1725  {
1726  //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1727 
1728  //fOutputAODBranch->SetOwner(kFALSE);
1729 
1730  AddAODBranch("TClonesArray", &fOutputAODBranch);
1731  }
1732 }
1733 
1734 //_______________________________________________________
1745 //_______________________________________________________
1747 {
1748  //-------
1749  // Step 1
1750 
1751  // Remove the contents of AOD branch output list set in the previous event
1752  fOutputAODBranch->Clear("C");
1753 
1754  LoadBranches();
1755 
1756  // Check if there is a centrality value, PbPb analysis, and if a centrality bin selection is requested
1757  // If we need a centrality bin, we select only those events in the corresponding bin.
1758  if( GetCentrality() && fCentralityBin[0] >= 0 && fCentralityBin[1] >= 0 )
1759  {
1760  Float_t cen = GetEventCentrality();
1761  if(cen > fCentralityBin[1] || cen < fCentralityBin[0]) return ; //reject events out of bin.
1762  }
1763 
1764  // Intermediate array with new clusters : init the array only once or clear from previous event
1765  if(!fCaloClusterArr) fCaloClusterArr = new TObjArray(10000);
1766  else fCaloClusterArr->Delete();//Clear("C"); it leaks?
1767 
1768  InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
1769 
1770  // Get the event, do some checks and settings
1771  CheckAndGetEvent() ;
1772 
1773  if (!fEvent)
1774  {
1775  AliDebug(1,Form("Skip Event %d", (Int_t) Entry()));
1776  return ;
1777  }
1778 
1779  // Init pointers, geometry, clusterizer, ocdb, aodb
1780 
1781  if(fAccessOCDB) AccessOCDB();
1782  if(fAccessOADB) AccessOADB(); // only once
1783 
1785 
1786  //-------
1787  // Step 2
1788 
1789  // Make clusters
1791  else ClusterizeCells() ;
1792 
1793  //-------
1794  // Step 3
1795 
1797 }
1798 
1799 
void FillAODCaloCells()
Put calo cells in standard branch.
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
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.
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
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.
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.
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.