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