AliPhysics  fe039ad (fe039ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // Properly instantiate the object
35 
40  AliEmcalContainer(),
41  fClusTimeCutLow(-10),
42  fClusTimeCutUp(10),
43  fExoticCut(kTRUE),
44  fDefaultClusterEnergy(-1),
45  fIncludePHOS(kFALSE),
46  fIncludePHOSonly(kFALSE),
47  fPhosMinNcells(0),
48  fPhosMinM02(0),
49  fEmcalMinM02(DBL_MIN),
50  fEmcalMaxM02(DBL_MAX),
51  fEmcalMaxM02CutEnergy(DBL_MAX)
52 {
53  fBaseClassName = "AliVCluster";
54  SetClassName("AliVCluster");
55 
56  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
57  fUserDefEnergyCut[i] = 0.;
58  }
59 }
60 
66  AliEmcalContainer(name),
67  fClusTimeCutLow(-10),
68  fClusTimeCutUp(10),
69  fExoticCut(kTRUE),
70  fDefaultClusterEnergy(-1),
71  fIncludePHOS(kFALSE),
72  fIncludePHOSonly(kFALSE),
73  fPhosMinNcells(0),
74  fPhosMinM02(0),
75  fEmcalMinM02(DBL_MIN),
76  fEmcalMaxM02(DBL_MAX),
77  fEmcalMaxM02CutEnergy(DBL_MAX)
78 {
79  fBaseClassName = "AliVCluster";
80  SetClassName("AliVCluster");
81 
82  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
83  fUserDefEnergyCut[i] = 0.;
84  }
85 }
86 
92 AliVCluster* AliClusterContainer::GetLeadingCluster(const char* opt)
93 {
94  TString option(opt);
95  option.ToLower();
96 
97  double (AliTLorentzVector::*momentum)() const = 0;
98 
99  if (option.Contains("e")) {
100  momentum = &AliTLorentzVector::E;
101  }
102  else {
103  momentum = &AliTLorentzVector::Et;
104  }
105 
107 
108  for (auto cluster : accepted_momentum()) {
109  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
110  clusterMax = cluster;
111  }
112  }
113 
114  return clusterMax.second;
115 }
116 
123 {
124  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
125  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
126  return vp;
127 
128 }
129 
136 {
137  AliVCluster *vc = GetCluster(i);
138  if (!vc) return 0;
139 
140  UInt_t rejectionReason = 0;
141  if (AcceptCluster(vc, rejectionReason))
142  return vc;
143  else {
144  AliDebug(2,"Cluster not accepted.");
145  return 0;
146  }
147 }
148 
155 {
156  Int_t i = GetIndexFromLabel(lab);
157  return GetCluster(i);
158 }
159 
166 {
167  Int_t i = GetIndexFromLabel(lab);
168  return GetAcceptCluster(i);
169 }
170 
176 {
177  const Int_t n = GetNEntries();
178  AliVCluster *c = 0;
179  do {
180  fCurrentID++;
181  if (fCurrentID >= n) break;
182  c = GetAcceptCluster(fCurrentID);
183  } while (!c);
184 
185  return c;
186 }
187 
193 {
194  const Int_t n = GetNEntries();
195  AliVCluster *c = 0;
196  do {
197  fCurrentID++;
198  if (fCurrentID >= n) break;
199  c = GetCluster(fCurrentID);
200  } while (!c);
201 
202  return c;
203 }
204 
205 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
206 {
207  if (mass < 0) mass = 0;
208 
209  Double_t energy = 0;
210 
211  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
212  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
213  }
214  else {
215  energy = vc->E();
216  }
217 
218  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
219 
220  Float_t pos[3];
221  vc->GetPosition(pos);
222 
223  pos[0]-=fVertex[0];
224  pos[1]-=fVertex[1];
225  pos[2]-=fVertex[2];
226 
227  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
228 
229  if (r > 1e-12) {
230  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
231  }
232  else {
233  AliInfo("Null cluster radius, momentum calculation not possible");
234  return kFALSE;
235  }
236 
237  return kTRUE;
238 }
239 
240 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
241 {
242  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
243 
244  if (vc) {
245  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
246  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
247  }
248  else {
249  vc->GetMomentum(mom, fVertex);
250  }
251  }
252  else {
253  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
254  return kFALSE;
255  }
256  return kTRUE;
257 }
258 
265 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
266 {
267  AliVCluster *vc = GetCluster(i);
268  return GetMomentum(mom, vc);
269 }
270 
277 {
278  AliVCluster *vc = GetNextCluster();
279  return GetMomentum(mom, vc);
280 }
281 
289 {
290  AliVCluster *vc = GetAcceptCluster(i);
291  return GetMomentum(mom, vc);
292 }
293 
300 {
301  AliVCluster *vc = GetNextAcceptCluster();
302  return GetMomentum(mom, vc);
303 }
304 
306 {
307  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
308  if (!r) return kFALSE;
309 
310  AliTLorentzVector mom;
311  GetMomentum(mom, i);
312 
313  return ApplyKinematicCuts(mom, rejectionReason);
314 }
315 
316 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
317 {
318  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
319  if (!r) return kFALSE;
320 
321  AliTLorentzVector mom;
322  GetMomentum(mom, clus);
323 
324  return ApplyKinematicCuts(mom, rejectionReason);
325 }
326 
332 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
333 {
334  if (!clus) {
335  rejectionReason |= kNullObject;
336  return kFALSE;
337  }
338 
339  Bool_t bInAcceptance = clus->IsEMCAL();
340  if (fIncludePHOS) {
341  bInAcceptance = clus->IsEMCAL() || (clus->GetType() == AliVCluster::kPHOSNeutral);
342  }
343  if (fIncludePHOSonly) {
344  bInAcceptance = (clus->GetType() == AliVCluster::kPHOSNeutral);
345  }
346  if (!bInAcceptance) {
347  rejectionReason |= kIsEMCalCut;
348  return kFALSE;
349  }
350 
351  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
352  rejectionReason |= kBitMapCut;
353  return kFALSE;
354  }
355 
356  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
357  rejectionReason |= kMCLabelCut;
358  return kFALSE;
359  }
360 
361  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
362  rejectionReason |= kMCLabelCut;
363  return kFALSE;
364  }
365 
366  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
367  rejectionReason |= kTimeCut;
368  return kFALSE;
369  }
370 
371  if (fExoticCut && clus->GetIsExotic()) {
372  rejectionReason |= kExoticCut;
373  return kFALSE;
374  }
375 
376  if (clus->IsEMCAL()) {
377  if (clus->GetNonLinCorrEnergy() < fEmcalMaxM02CutEnergy) {
378  if(clus->GetM02() < fEmcalMinM02 || clus->GetM02() > fEmcalMaxM02) {
379  rejectionReason |= kExoticCut; // Not really true, but there is a lack of leftover bits
380  return kFALSE;
381  }
382  }
383  }
384 
385  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
386  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
387  rejectionReason |= kEnergyCut;
388  return kFALSE;
389  }
390  }
391 
393  if (clus->GetType() == AliVCluster::kPHOSNeutral) {
394  if (clus->GetNCells() < fPhosMinNcells) {
395  rejectionReason |= kExoticCut;
396  return kFALSE;
397  }
398 
399  if (clus->GetM02() < fPhosMinM02) {
400  rejectionReason |= kExoticCut;
401  return kFALSE;
402  }
403  }
404  }
405 
406  return kTRUE;
407 }
408 
414 {
415  UInt_t rejectionReason = 0;
416  Int_t nClus = 0;
417  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
418  AliVCluster *clust = this->GetCluster(iclust);
419  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
420  }
421  return nClus;
422 }
423 
430 {
431  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
432  return fUserDefEnergyCut[t];
433  }
434  else {
435  return fMinE;
436  }
437 }
438 
445 {
446  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
447  fUserDefEnergyCut[t] = cut;
448  }
449  else {
450  fMinE = cut;
451  }
452 }
453 
462 void AliClusterContainer::SetArray(const AliVEvent * event)
463 {
464  AliEmcalContainer::SetArray(event);
465 
466  // Register TClonesArray in index map
468 }
469 
476  return AliClusterIterableContainer(this, false);
477 }
478 
485  return AliClusterIterableContainer(this, true);
486 }
487 
494  return AliClusterIterableMomentumContainer(this, false);
495 }
496 
503  return AliClusterIterableMomentumContainer(this, true);
504 }
505 
506 const char* AliClusterContainer::GetTitle() const
507 {
508  static TString clusterString;
509 
511 
512  if (Ecut == 0) {
513  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
514  }
515  else if (Ecut < 1.0) {
516  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
517  }
518  else {
519  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
520  }
521 
522  return clusterString.Data();
523 }
524 
525 TString AliClusterContainer::GetDefaultArrayName(const AliVEvent * const ev) const {
526  if(ev->IsA() == AliAODEvent::Class()) return "caloClusters";
527  else if(ev->IsA() == AliESDEvent::Class()) return "CaloClusters";
528  else return "";
529 }
530 
531 
532 /******************************************
533  * Unit tests *
534  ******************************************/
535 
536 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
537  std::vector<AliVCluster *> reference, variation;
538  AliVCluster *test = NULL;
539  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
540  test = cont->GetCluster(iclust);
541  if(!iteratorType){
542  UInt_t rejectionReason = 0;
543  if(!cont->AcceptCluster(test, rejectionReason)) continue;
544  }
545  reference.push_back(test);
546  }
547 
548  if(!iteratorType){
549  // test accept iterator
550  for(auto cluster : cont->accepted()){
551  variation.push_back(cluster);
552  }
553  } else {
554  // test all iterator
555  for(auto cluster : cont->all()){
556  variation.push_back(cluster);
557  }
558  }
559 
560  if(verbose){
561  // Printing cluster addresses:
562  if(reference.size() < 30){
563  std::cout << "Clusters in reference container: " << std::endl;
564  std::cout << "===========================================" << std::endl;
565  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
566  std::cout << "Address: " << *refit << std::endl;
567  }
568  }
569  if(variation.size() < 30){
570  std::cout << "Clusters in test container: " << std::endl;
571  std::cout << "===========================================" << std::endl;
572  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
573  std::cout << "Address: " << *varit << std::endl;
574  }
575  }
576  }
577 
578  int testresult = 0;
579  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
580  bool failure = false;
581  // first test: cleck if all clusters are found by the iterator
582  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
583  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
584  if(verbose)
585  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
586  failure = true;
587  break;
588  }
589  }
590  if(!failure){
591  // second test: check if there are no extra clusters found
592  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
593  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
594  if(verbose)
595  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
596  failure = true;
597  break;
598  }
599  }
600  if(failure) testresult = 2;
601  } else {
602  testresult = 1;
603  }
604 
605  if(verbose){
606  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
607  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
608  std::cout << "Test result: " << testresult << std::endl;
609  std::cout << "-----------------------------------------------------------------------" << std::endl;
610  }
611  return testresult;
612 }
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.
AliVCluster * GetNextAcceptCluster()
Int_t fPhosMinNcells
min number of phos cells per cluster