AliPhysics  master (3d17d9d)
AliClusterContainer.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 #include <algorithm>
16 #include <cfloat>
17 #include <iostream>
18 #include <vector>
19 #include <TClonesArray.h>
20 
21 #include "AliAODEvent.h"
22 #include "AliESDEvent.h"
23 #include "AliVCaloCells.h"
24 #include "AliVEvent.h"
25 #include "AliLog.h"
26 #include "AliTLorentzVector.h"
27 
28 #include "AliClusterContainer.h"
29 
31 ClassImp(AliClusterContainer);
33 
34 // string to enum map for use with the %YAML config
35 const std::map <std::string, AliVCluster::VCluUserDefEnergy_t> AliClusterContainer::fgkClusterEnergyTypeMap = {
36  {"kNonLinCorr", AliVCluster::kNonLinCorr },
37  {"kHadCorr", AliVCluster::kHadCorr },
38  {"kUserDefEnergy1", AliVCluster::kUserDefEnergy1 },
39  {"kUserDefEnergy2", AliVCluster::kUserDefEnergy2 }
40 };
41 
42 // Properly instantiate the object
44 
50  fEMCALCells(nullptr),
51  fClusTimeCutLow(-10),
52  fClusTimeCutUp(10),
53  fExoticCut(kTRUE),
54  fDefaultClusterEnergy(-1),
55  fIncludePHOS(kFALSE),
56  fIncludePHOSonly(kFALSE),
57  fPhosMinNcells(0),
58  fPhosMinM02(0),
59  fEmcalMinM02(-1.),
60  fEmcalMaxM02(DBL_MAX),
61  fEmcalMaxM02CutEnergy(DBL_MAX),
62  fMaxFracEnergyLeadingCell(DBL_MAX)
63 {
64  fBaseClassName = "AliVCluster";
65  SetClassName("AliVCluster");
66 
67  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
68  fUserDefEnergyCut[i] = 0.;
69  }
70 }
71 
77  AliEmcalContainer(name),
79  fClusTimeCutLow(-10),
80  fClusTimeCutUp(10),
81  fExoticCut(kTRUE),
83  fIncludePHOS(kFALSE),
84  fIncludePHOSonly(kFALSE),
85  fPhosMinNcells(0),
86  fPhosMinM02(0),
87  fEmcalMinM02(-1.),
88  fEmcalMaxM02(DBL_MAX),
89  fEmcalMaxM02CutEnergy(DBL_MAX),
91 {
92  fBaseClassName = "AliVCluster";
93  SetClassName("AliVCluster");
94 
95  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
96  fUserDefEnergyCut[i] = 0.;
97  }
98 }
99 
106 {
107  TString option(opt);
108  option.ToLower();
109 
110  double (AliTLorentzVector::*momentum)() const = 0;
111 
112  if (option.Contains("e")) {
113  momentum = &AliTLorentzVector::E;
114  }
115  else {
116  momentum = &AliTLorentzVector::Et;
117  }
118 
120 
121  for (auto cluster : accepted_momentum()) {
122  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
123  clusterMax = cluster;
124  }
125  }
126 
127  return clusterMax.second;
128 }
129 
136 {
137  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
138  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
139  return vp;
140 
141 }
142 
149 {
150  AliVCluster *vc = GetCluster(i);
151  if (!vc) return 0;
152 
153  UInt_t rejectionReason = 0;
154  if (AcceptCluster(vc, rejectionReason))
155  return vc;
156  else {
157  AliDebug(2,"Cluster not accepted.");
158  return 0;
159  }
160 }
161 
168 {
169  Int_t i = GetIndexFromLabel(lab);
170  return GetCluster(i);
171 }
172 
179 {
180  Int_t i = GetIndexFromLabel(lab);
181  return GetAcceptCluster(i);
182 }
183 
189 {
190  const Int_t n = GetNEntries();
191  AliVCluster *c = 0;
192  do {
193  fCurrentID++;
194  if (fCurrentID >= n) break;
196  } while (!c);
197 
198  return c;
199 }
200 
206 {
207  const Int_t n = GetNEntries();
208  AliVCluster *c = 0;
209  do {
210  fCurrentID++;
211  if (fCurrentID >= n) break;
212  c = GetCluster(fCurrentID);
213  } while (!c);
214 
215  return c;
216 }
217 
218 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
219 {
220  if (mass < 0) mass = 0;
221 
222  Double_t energy = 0;
223 
224  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
225  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
226  }
227  else {
228  energy = vc->E();
229  }
230 
231  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
232 
233  Float_t pos[3];
234  vc->GetPosition(pos);
235 
236  pos[0]-=fVertex[0];
237  pos[1]-=fVertex[1];
238  pos[2]-=fVertex[2];
239 
240  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
241 
242  if (r > 1e-12) {
243  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
244  }
245  else {
246  AliInfo("Null cluster radius, momentum calculation not possible");
247  return kFALSE;
248  }
249 
250  return kTRUE;
251 }
252 
253 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
254 {
255  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
256 
257  if (vc) {
258  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
259  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
260  }
261  else {
262  vc->GetMomentum(mom, fVertex);
263  }
264  }
265  else {
266  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
267  return kFALSE;
268  }
269  return kTRUE;
270 }
271 
278 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
279 {
280  AliVCluster *vc = GetCluster(i);
281  return GetMomentum(mom, vc);
282 }
283 
290 {
291  AliVCluster *vc = GetNextCluster();
292  return GetMomentum(mom, vc);
293 }
294 
302 {
303  AliVCluster *vc = GetAcceptCluster(i);
304  return GetMomentum(mom, vc);
305 }
306 
313 {
314  AliVCluster *vc = GetNextAcceptCluster();
315  return GetMomentum(mom, vc);
316 }
317 
319 {
320  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
321  if (!r) return kFALSE;
322 
323  AliTLorentzVector mom;
324  GetMomentum(mom, i);
325 
326  return ApplyKinematicCuts(mom, rejectionReason);
327 }
328 
329 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
330 {
331  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
332  if (!r) return kFALSE;
333 
334  AliTLorentzVector mom;
335  GetMomentum(mom, clus);
336 
337  return ApplyKinematicCuts(mom, rejectionReason);
338 }
339 
346 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
347 {
348  if (!clus) {
349  rejectionReason |= kNullObject;
350  return kFALSE;
351  }
352 
353  Bool_t bInAcceptance = clus->IsEMCAL();
354  if (fIncludePHOS) {
355  bInAcceptance = clus->IsEMCAL() || (clus->GetType() == AliVCluster::kPHOSNeutral);
356  }
357  if (fIncludePHOSonly) {
358  bInAcceptance = (clus->GetType() == AliVCluster::kPHOSNeutral);
359  }
360  if (!bInAcceptance) {
361  rejectionReason |= kIsEMCalCut;
362  return kFALSE;
363  }
364 
365  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
366  rejectionReason |= kBitMapCut;
367  return kFALSE;
368  }
369 
370  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
371  rejectionReason |= kMCLabelCut;
372  return kFALSE;
373  }
374 
375  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
376  rejectionReason |= kMCLabelCut;
377  return kFALSE;
378  }
379 
380  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
381  rejectionReason |= kTimeCut;
382  return kFALSE;
383  }
384 
385  if (fExoticCut && clus->GetIsExotic()) {
386  rejectionReason |= kExoticCut;
387  return kFALSE;
388  }
389 
390  if (clus->IsEMCAL()) {
391  if (clus->GetNonLinCorrEnergy() < fEmcalMaxM02CutEnergy) {
392  if(clus->GetM02() < fEmcalMinM02 || clus->GetM02() > fEmcalMaxM02) {
393  rejectionReason |= kExoticCut; // Not really true, but there is a lack of leftover bits
394  return kFALSE;
395  }
396  }
397 
398  }
399 
400  if(fEMCALCells) {
401  double ecellmax = 0;
402  for(int icell = 0; icell < clus->GetNCells(); icell++){
403  double celltmp = fEMCALCells->GetCellAmplitude(clus->GetCellAbsId(icell));
404  if(celltmp > ecellmax) ecellmax = celltmp;
405  }
406  if(ecellmax/clus->E() > fMaxFracEnergyLeadingCell) {
407  rejectionReason |= kExoticCut;
408  return kFALSE;
409  }
410  }
411 
412  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
413  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
414  rejectionReason |= kEnergyCut;
415  return kFALSE;
416  }
417  }
418 
420  if (clus->GetType() == AliVCluster::kPHOSNeutral) {
421  if (clus->GetNCells() < fPhosMinNcells) {
422  rejectionReason |= kExoticCut;
423  return kFALSE;
424  }
425 
426  if (clus->GetM02() < fPhosMinM02) {
427  rejectionReason |= kExoticCut;
428  return kFALSE;
429  }
430  }
431  }
432 
433  return kTRUE;
434 }
435 
441 {
442  UInt_t rejectionReason = 0;
443  Int_t nClus = 0;
444  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
445  AliVCluster *clust = this->GetCluster(iclust);
446  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
447  }
448  return nClus;
449 }
450 
457 {
458  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
459  return fUserDefEnergyCut[t];
460  }
461  else {
462  return fMinE;
463  }
464 }
465 
472 {
473  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
474  fUserDefEnergyCut[t] = cut;
475  }
476  else {
477  fMinE = cut;
478  }
479 }
480 
489 void AliClusterContainer::SetArray(const AliVEvent * event)
490 {
492 
493  // Connect EMCAL cells
494  fEMCALCells = event->GetEMCALCells();
495 
496  // Register TClonesArray in index map
498 }
499 
506  return AliClusterIterableContainer(this, false);
507 }
508 
515  return AliClusterIterableContainer(this, true);
516 }
517 
524  return AliClusterIterableMomentumContainer(this, false);
525 }
526 
533  return AliClusterIterableMomentumContainer(this, true);
534 }
535 
536 const char* AliClusterContainer::GetTitle() const
537 {
538  static TString clusterString;
539 
541 
542  if (Ecut == 0) {
543  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
544  }
545  else if (Ecut < 1.0) {
546  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
547  }
548  else {
549  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
550  }
551 
552  return clusterString.Data();
553 }
554 
555 TString AliClusterContainer::GetDefaultArrayName(const AliVEvent * const ev) const {
556  if(ev->IsA() == AliAODEvent::Class()) return "caloClusters";
557  else if(ev->IsA() == AliESDEvent::Class()) return "CaloClusters";
558  else return "";
559 }
560 
561 
562 /******************************************
563  * Unit tests *
564  ******************************************/
565 
566 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
567  std::vector<AliVCluster *> reference, variation;
568  AliVCluster *test = NULL;
569  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
570  test = cont->GetCluster(iclust);
571  if(!iteratorType){
572  UInt_t rejectionReason = 0;
573  if(!cont->AcceptCluster(test, rejectionReason)) continue;
574  }
575  reference.push_back(test);
576  }
577 
578  if(!iteratorType){
579  // test accept iterator
580  for(auto cluster : cont->accepted()){
581  variation.push_back(cluster);
582  }
583  } else {
584  // test all iterator
585  for(auto cluster : cont->all()){
586  variation.push_back(cluster);
587  }
588  }
589 
590  if(verbose){
591  // Printing cluster addresses:
592  if(reference.size() < 30){
593  std::cout << "Clusters in reference container: " << std::endl;
594  std::cout << "===========================================" << std::endl;
595  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
596  std::cout << "Address: " << *refit << std::endl;
597  }
598  }
599  if(variation.size() < 30){
600  std::cout << "Clusters in test container: " << std::endl;
601  std::cout << "===========================================" << std::endl;
602  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
603  std::cout << "Address: " << *varit << std::endl;
604  }
605  }
606  }
607 
608  int testresult = 0;
609  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
610  bool failure = false;
611  // first test: cleck if all clusters are found by the iterator
612  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
613  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
614  if(verbose)
615  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
616  failure = true;
617  break;
618  }
619  }
620  if(!failure){
621  // second test: check if there are no extra clusters found
622  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
623  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
624  if(verbose)
625  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
626  failure = true;
627  break;
628  }
629  }
630  if(failure) testresult = 2;
631  } else {
632  testresult = 1;
633  }
634 
635  if(verbose){
636  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
637  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
638  std::cout << "Test result: " << testresult << std::endl;
639  std::cout << "-----------------------------------------------------------------------" << std::endl;
640  }
641  return testresult;
642 }
AliVCluster * GetAcceptClusterWithLabel(Int_t lab) const
virtual TString GetDefaultArrayName(const AliVEvent *const ev) const
virtual Bool_t ApplyClusterCuts(const AliVCluster *clus, UInt_t &rejectionReason) const
Bool_t GetNextAcceptMomentum(TLorentzVector &mom)
virtual Bool_t ApplyKinematicCuts(const AliTLorentzVector &mom, UInt_t &rejectionReason) const
Apply kinematical selection to the momentum vector provided.
const char * GetTitle() const
double Double_t
Definition: External.C:58
Double_t fMinE
Max. cut on particle energy.
EMCALIterableContainer::AliEmcalIterableContainerT< AliVCluster, EMCALIterableContainer::operator_star_object< AliVCluster > > AliClusterIterableContainer
Bool_t GetNextMomentum(TLorentzVector &mom)
Int_t GetIndexFromLabel(Int_t lab) const
Get the index in the container from a given label.
const AliClusterIterableMomentumContainer accepted_momentum() const
Double_t mass
energy
Definition: HFPtSpectrum.C:45
Bool_t fIncludePHOS
flag to accept PHOS clusters in addition to EMCal clusters
Declaration of class AliTLorentzVector.
UInt_t fBitMap
bitmap mask
Double_t fClusTimeCutLow
low time cut for clusters
Double_t fUserDefEnergyCut[AliVCluster::kLastUserDefEnergy+1]
cut on the energy of the cluster after higher level corrections (see AliVCluster.h) ...
TCanvas * c
Definition: TestFitELoss.C:172
void SetClusUserDefEnergyCut(Int_t t, Double_t cut)
AliVCluster * GetLeadingCluster(const char *opt="")
Cell time cut not passed.
Double_t GetClusUserDefEnergyCut(Int_t t) const
Int_t GetDefaultClusterEnergy() const
Double_t fMaxFracEnergyLeadingCell
max fraction of energy in the leading cell
const AliClusterIterableContainer all() const
enum AliVCluster::VCluUserDefEnergy_t VCluUserDefEnergy_t
Cluster is exotic cluster.
Double_t fEmcalMaxM02CutEnergy
max EMCal cluster energy for which to apply M02 cut
std::pair< AliTLorentzVector, T * > momentum_object_pair
Energy below threshold.
int Int_t
Definition: External.C:63
Double_t fEmcalMaxM02
max value of M02 for EMCAL clusters
Double_t fEmcalMinM02
min value of M02 for EMCAL clusters
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Base class for container structures within the EMCAL framework.
virtual Bool_t AcceptCluster(Int_t i, UInt_t &rejectionReason) const
AliVCluster * GetAcceptCluster(Int_t i) const
Int_t GetNClusters() const
Int_t GetNAcceptedClusters() const
const AliClusterIterableContainer accepted() const
TClonesArray * GetArray() const
Double_t fClusTimeCutUp
up time cut for clusters
AliVCluster * GetCluster(Int_t i) const
void SetClassName(const char *clname)
AliVCaloCells * fEMCALCells
pointer to EMCAL cells object
Bool_t fExoticCut
reject clusters marked as "exotic"
AliVCluster * GetNextCluster()
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Int_t fMinMCLabel
minimum MC label
EMCALIterableContainer::AliEmcalIterableContainerT< AliVCluster, EMCALIterableContainer::operator_star_pair< AliVCluster > > AliClusterIterableMomentumContainer
const TString & GetArrayName() const
Double_t fMassHypothesis
if < 0 it will use a PID mass when available
Cluster not in the EMCAL.
const AliClusterIterableMomentumContainer all_momentum() const
Double_t fVertex[3]
! event vertex array
Bool_t GetAcceptMomentum(TLorentzVector &mom, Int_t i) const
AliVCluster * GetClusterWithLabel(Int_t lab) const
Double_t fPhosMinM02
min value of M02 for phos clusters
int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose)
virtual void SetArray(const AliVEvent *event)
TString fBaseClassName
name of the base class that this container can handle
Bool_t GetMomentum(TLorentzVector &mom, const AliVCluster *vc, Double_t mass) const
Int_t fMaxMCLabel
maximum MC label
void SetArray(const AliVEvent *event)
void test(int runnumber=195345)
Bool_t fIncludePHOSonly
flag to accept only PHOS clusters (and reject EMCal clusters)
TClonesArray * fClArray
! Pointer to array in input event
Int_t GetNEntries() const
bool Bool_t
Definition: External.C:53
Int_t fDefaultClusterEnergy
default cluster energy: -1 for clus->E(); otherwise clus->GetUserDefEnergy(fDefaultClusterEnergy) ...
Int_t fCurrentID
! current ID for automatic loops
static AliEmcalContainerIndexMap< TClonesArray, AliVCluster > fgEmcalContainerIndexMap
! Mapping from containers to indices
Container structure for EMCAL clusters.
static const std::map< std::string, VCluUserDefEnergy_t > fgkClusterEnergyTypeMap
Relates string to the cluster energy enumeration for YAML configuration.
AliVCluster * GetNextAcceptCluster()
Int_t fPhosMinNcells
min number of phos cells per cluster