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