AliPhysics  aaf9c62 (aaf9c62)
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 "AliVEvent.h"
24 #include "AliLog.h"
25 #include "AliTLorentzVector.h"
26 
27 #include "AliClusterContainer.h"
28 
30 ClassImp(AliClusterContainer);
32 
33 // string to enum map for use with the %YAML config
34 const std::map <std::string, AliVCluster::VCluUserDefEnergy_t> AliClusterContainer::fgkClusterEnergyTypeMap = {
35  {"kNonLinCorr", AliVCluster::kNonLinCorr },
36  {"kHadCorr", AliVCluster::kHadCorr },
37  {"kUserDefEnergy1", AliVCluster::kUserDefEnergy1 },
38  {"kUserDefEnergy2", AliVCluster::kUserDefEnergy2 }
39 };
40 
41 // Properly instantiate the object
43 
48  AliEmcalContainer(),
49  fClusTimeCutLow(-10),
50  fClusTimeCutUp(10),
51  fExoticCut(kTRUE),
52  fDefaultClusterEnergy(-1),
53  fIncludePHOS(kFALSE),
54  fIncludePHOSonly(kFALSE),
55  fPhosMinNcells(0),
56  fPhosMinM02(0),
57  fEmcalMinM02(-1.),
58  fEmcalMaxM02(DBL_MAX),
59  fEmcalMaxM02CutEnergy(DBL_MAX)
60 {
61  fBaseClassName = "AliVCluster";
62  SetClassName("AliVCluster");
63 
64  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
65  fUserDefEnergyCut[i] = 0.;
66  }
67 }
68 
74  AliEmcalContainer(name),
75  fClusTimeCutLow(-10),
76  fClusTimeCutUp(10),
77  fExoticCut(kTRUE),
79  fIncludePHOS(kFALSE),
80  fIncludePHOSonly(kFALSE),
81  fPhosMinNcells(0),
82  fPhosMinM02(0),
83  fEmcalMinM02(-1.),
84  fEmcalMaxM02(DBL_MAX),
85  fEmcalMaxM02CutEnergy(DBL_MAX)
86 {
87  fBaseClassName = "AliVCluster";
88  SetClassName("AliVCluster");
89 
90  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
91  fUserDefEnergyCut[i] = 0.;
92  }
93 }
94 
101 {
102  TString option(opt);
103  option.ToLower();
104 
105  double (AliTLorentzVector::*momentum)() const = 0;
106 
107  if (option.Contains("e")) {
108  momentum = &AliTLorentzVector::E;
109  }
110  else {
111  momentum = &AliTLorentzVector::Et;
112  }
113 
115 
116  for (auto cluster : accepted_momentum()) {
117  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
118  clusterMax = cluster;
119  }
120  }
121 
122  return clusterMax.second;
123 }
124 
131 {
132  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
133  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
134  return vp;
135 
136 }
137 
144 {
145  AliVCluster *vc = GetCluster(i);
146  if (!vc) return 0;
147 
148  UInt_t rejectionReason = 0;
149  if (AcceptCluster(vc, rejectionReason))
150  return vc;
151  else {
152  AliDebug(2,"Cluster not accepted.");
153  return 0;
154  }
155 }
156 
163 {
164  Int_t i = GetIndexFromLabel(lab);
165  return GetCluster(i);
166 }
167 
174 {
175  Int_t i = GetIndexFromLabel(lab);
176  return GetAcceptCluster(i);
177 }
178 
184 {
185  const Int_t n = GetNEntries();
186  AliVCluster *c = 0;
187  do {
188  fCurrentID++;
189  if (fCurrentID >= n) break;
190  c = GetAcceptCluster(fCurrentID);
191  } while (!c);
192 
193  return c;
194 }
195 
201 {
202  const Int_t n = GetNEntries();
203  AliVCluster *c = 0;
204  do {
205  fCurrentID++;
206  if (fCurrentID >= n) break;
207  c = GetCluster(fCurrentID);
208  } while (!c);
209 
210  return c;
211 }
212 
213 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
214 {
215  if (mass < 0) mass = 0;
216 
217  Double_t energy = 0;
218 
219  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
220  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
221  }
222  else {
223  energy = vc->E();
224  }
225 
226  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
227 
228  Float_t pos[3];
229  vc->GetPosition(pos);
230 
231  pos[0]-=fVertex[0];
232  pos[1]-=fVertex[1];
233  pos[2]-=fVertex[2];
234 
235  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
236 
237  if (r > 1e-12) {
238  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
239  }
240  else {
241  AliInfo("Null cluster radius, momentum calculation not possible");
242  return kFALSE;
243  }
244 
245  return kTRUE;
246 }
247 
248 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
249 {
250  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
251 
252  if (vc) {
253  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
254  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
255  }
256  else {
257  vc->GetMomentum(mom, fVertex);
258  }
259  }
260  else {
261  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
262  return kFALSE;
263  }
264  return kTRUE;
265 }
266 
273 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
274 {
275  AliVCluster *vc = GetCluster(i);
276  return GetMomentum(mom, vc);
277 }
278 
285 {
286  AliVCluster *vc = GetNextCluster();
287  return GetMomentum(mom, vc);
288 }
289 
297 {
298  AliVCluster *vc = GetAcceptCluster(i);
299  return GetMomentum(mom, vc);
300 }
301 
308 {
309  AliVCluster *vc = GetNextAcceptCluster();
310  return GetMomentum(mom, vc);
311 }
312 
314 {
315  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
316  if (!r) return kFALSE;
317 
318  AliTLorentzVector mom;
319  GetMomentum(mom, i);
320 
321  return ApplyKinematicCuts(mom, rejectionReason);
322 }
323 
324 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
325 {
326  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
327  if (!r) return kFALSE;
328 
329  AliTLorentzVector mom;
330  GetMomentum(mom, clus);
331 
332  return ApplyKinematicCuts(mom, rejectionReason);
333 }
334 
341 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
342 {
343  if (!clus) {
344  rejectionReason |= kNullObject;
345  return kFALSE;
346  }
347 
348  Bool_t bInAcceptance = clus->IsEMCAL();
349  if (fIncludePHOS) {
350  bInAcceptance = clus->IsEMCAL() || (clus->GetType() == AliVCluster::kPHOSNeutral);
351  }
352  if (fIncludePHOSonly) {
353  bInAcceptance = (clus->GetType() == AliVCluster::kPHOSNeutral);
354  }
355  if (!bInAcceptance) {
356  rejectionReason |= kIsEMCalCut;
357  return kFALSE;
358  }
359 
360  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
361  rejectionReason |= kBitMapCut;
362  return kFALSE;
363  }
364 
365  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
366  rejectionReason |= kMCLabelCut;
367  return kFALSE;
368  }
369 
370  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
371  rejectionReason |= kMCLabelCut;
372  return kFALSE;
373  }
374 
375  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
376  rejectionReason |= kTimeCut;
377  return kFALSE;
378  }
379 
380  if (fExoticCut && clus->GetIsExotic()) {
381  rejectionReason |= kExoticCut;
382  return kFALSE;
383  }
384 
385  if (clus->IsEMCAL()) {
386  if (clus->GetNonLinCorrEnergy() < fEmcalMaxM02CutEnergy) {
387  if(clus->GetM02() < fEmcalMinM02 || clus->GetM02() > fEmcalMaxM02) {
388  rejectionReason |= kExoticCut; // Not really true, but there is a lack of leftover bits
389  return kFALSE;
390  }
391  }
392  }
393 
394  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
395  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
396  rejectionReason |= kEnergyCut;
397  return kFALSE;
398  }
399  }
400 
402  if (clus->GetType() == AliVCluster::kPHOSNeutral) {
403  if (clus->GetNCells() < fPhosMinNcells) {
404  rejectionReason |= kExoticCut;
405  return kFALSE;
406  }
407 
408  if (clus->GetM02() < fPhosMinM02) {
409  rejectionReason |= kExoticCut;
410  return kFALSE;
411  }
412  }
413  }
414 
415  return kTRUE;
416 }
417 
423 {
424  UInt_t rejectionReason = 0;
425  Int_t nClus = 0;
426  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
427  AliVCluster *clust = this->GetCluster(iclust);
428  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
429  }
430  return nClus;
431 }
432 
439 {
440  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
441  return fUserDefEnergyCut[t];
442  }
443  else {
444  return fMinE;
445  }
446 }
447 
454 {
455  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
456  fUserDefEnergyCut[t] = cut;
457  }
458  else {
459  fMinE = cut;
460  }
461 }
462 
471 void AliClusterContainer::SetArray(const AliVEvent * event)
472 {
473  AliEmcalContainer::SetArray(event);
474 
475  // Register TClonesArray in index map
477 }
478 
485  return AliClusterIterableContainer(this, false);
486 }
487 
494  return AliClusterIterableContainer(this, true);
495 }
496 
503  return AliClusterIterableMomentumContainer(this, false);
504 }
505 
512  return AliClusterIterableMomentumContainer(this, true);
513 }
514 
515 const char* AliClusterContainer::GetTitle() const
516 {
517  static TString clusterString;
518 
520 
521  if (Ecut == 0) {
522  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
523  }
524  else if (Ecut < 1.0) {
525  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
526  }
527  else {
528  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
529  }
530 
531  return clusterString.Data();
532 }
533 
534 TString AliClusterContainer::GetDefaultArrayName(const AliVEvent * const ev) const {
535  if(ev->IsA() == AliAODEvent::Class()) return "caloClusters";
536  else if(ev->IsA() == AliESDEvent::Class()) return "CaloClusters";
537  else return "";
538 }
539 
540 
541 /******************************************
542  * Unit tests *
543  ******************************************/
544 
545 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
546  std::vector<AliVCluster *> reference, variation;
547  AliVCluster *test = NULL;
548  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
549  test = cont->GetCluster(iclust);
550  if(!iteratorType){
551  UInt_t rejectionReason = 0;
552  if(!cont->AcceptCluster(test, rejectionReason)) continue;
553  }
554  reference.push_back(test);
555  }
556 
557  if(!iteratorType){
558  // test accept iterator
559  for(auto cluster : cont->accepted()){
560  variation.push_back(cluster);
561  }
562  } else {
563  // test all iterator
564  for(auto cluster : cont->all()){
565  variation.push_back(cluster);
566  }
567  }
568 
569  if(verbose){
570  // Printing cluster addresses:
571  if(reference.size() < 30){
572  std::cout << "Clusters in reference container: " << std::endl;
573  std::cout << "===========================================" << std::endl;
574  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
575  std::cout << "Address: " << *refit << std::endl;
576  }
577  }
578  if(variation.size() < 30){
579  std::cout << "Clusters in test container: " << std::endl;
580  std::cout << "===========================================" << std::endl;
581  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
582  std::cout << "Address: " << *varit << std::endl;
583  }
584  }
585  }
586 
587  int testresult = 0;
588  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
589  bool failure = false;
590  // first test: cleck if all clusters are found by the iterator
591  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
592  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
593  if(verbose)
594  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
595  failure = true;
596  break;
597  }
598  }
599  if(!failure){
600  // second test: check if there are no extra clusters found
601  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
602  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
603  if(verbose)
604  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
605  failure = true;
606  break;
607  }
608  }
609  if(failure) testresult = 2;
610  } else {
611  testresult = 1;
612  }
613 
614  if(verbose){
615  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
616  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
617  std::cout << "Test result: " << testresult << std::endl;
618  std::cout << "-----------------------------------------------------------------------" << std::endl;
619  }
620  return testresult;
621 }
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)
const char * GetTitle() const
double Double_t
Definition: External.C:58
EMCALIterableContainer::AliEmcalIterableContainerT< AliVCluster, EMCALIterableContainer::operator_star_object< AliVCluster > > AliClusterIterableContainer
Bool_t GetNextMomentum(TLorentzVector &mom)
const AliClusterIterableMomentumContainer accepted_momentum() const
Double_t mass
energy
Definition: HFPtSpectrum.C:44
Bool_t fIncludePHOS
flag to accept PHOS clusters in addition to EMCal clusters
Declaration of class AliTLorentzVector.
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="")
Double_t GetClusUserDefEnergyCut(Int_t t) const
Int_t GetDefaultClusterEnergy() const
const AliClusterIterableContainer all() const
enum AliVCluster::VCluUserDefEnergy_t VCluUserDefEnergy_t
Double_t fEmcalMaxM02CutEnergy
max EMCal cluster energy for which to apply M02 cut
std::pair< AliTLorentzVector, T * > momentum_object_pair
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
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
Double_t fClusTimeCutUp
up time cut for clusters
AliVCluster * GetCluster(Int_t i) const
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)
EMCALIterableContainer::AliEmcalIterableContainerT< AliVCluster, EMCALIterableContainer::operator_star_pair< AliVCluster > > AliClusterIterableMomentumContainer
const AliClusterIterableMomentumContainer all_momentum() const
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)
Bool_t GetMomentum(TLorentzVector &mom, const AliVCluster *vc, Double_t mass) const
void SetArray(const AliVEvent *event)
void test(int runnumber=195345)
Bool_t fIncludePHOSonly
flag to accept only PHOS clusters (and reject EMCal clusters)
bool Bool_t
Definition: External.C:53
Int_t fDefaultClusterEnergy
default cluster energy: -1 for clus->E(); otherwise clus->GetUserDefEnergy(fDefaultClusterEnergy) ...
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