AliPhysics  8bb951a (8bb951a)
 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 <iostream>
17 #include <vector>
18 #include <TClonesArray.h>
19 
20 #include "AliVEvent.h"
21 #include "AliLog.h"
22 #include "AliTLorentzVector.h"
23 
24 #include "AliClusterContainer.h"
25 
29 
35  fClusTimeCutLow(-10),
36  fClusTimeCutUp(10),
37  fExoticCut(kTRUE),
38  fDefaultClusterEnergy(-1)
39 {
40  fBaseClassName = "AliVCluster";
41  SetClassName("AliVCluster");
42 
43  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
44  fUserDefEnergyCut[i] = 0.;
45  }
46 }
47 
53  AliEmcalContainer(name),
54  fClusTimeCutLow(-10),
55  fClusTimeCutUp(10),
56  fExoticCut(kTRUE),
57  fDefaultClusterEnergy(-1)
58 {
59  fBaseClassName = "AliVCluster";
60  SetClassName("AliVCluster");
61 
62  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
63  fUserDefEnergyCut[i] = 0.;
64  }
65 }
66 
72 AliVCluster* AliClusterContainer::GetLeadingCluster(const char* opt)
73 {
74  TString option(opt);
75  option.ToLower();
76 
77  Int_t tempID = fCurrentID;
79 
80  AliVCluster *clusterMax = GetNextAcceptCluster();
81  AliVCluster *cluster = 0;
82 
83  if (option.Contains("e")) {
84  while ((cluster = GetNextAcceptCluster())) {
85  if (cluster->E() > clusterMax->E()) clusterMax = cluster;
86  }
87  }
88  else {
89  Double_t et = 0;
90  Double_t etmax = 0;
91  while ((cluster = GetNextAcceptCluster())) {
92  TLorentzVector mom;
93  cluster->GetMomentum(mom,const_cast<Double_t*>(fVertex));
94  et = mom.Et();
95  if (et > etmax) {
96  clusterMax = cluster;
97  etmax = et;
98  }
99  }
100  }
101 
102  fCurrentID = tempID;
103 
104  return clusterMax;
105 }
106 
112 AliVCluster* AliClusterContainer::GetCluster(Int_t i) const
113 {
114  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
115  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
116  return vp;
117 
118 }
119 
125 AliVCluster* AliClusterContainer::GetAcceptCluster(Int_t i) const
126 {
127  AliVCluster *vc = GetCluster(i);
128  if (!vc) return 0;
129 
130  UInt_t rejectionReason = 0;
131  if (AcceptCluster(vc, rejectionReason))
132  return vc;
133  else {
134  AliDebug(2,"Cluster not accepted.");
135  return 0;
136  }
137 }
138 
144 AliVCluster* AliClusterContainer::GetClusterWithLabel(Int_t lab) const
145 {
146  Int_t i = GetIndexFromLabel(lab);
147  return GetCluster(i);
148 }
149 
156 {
157  Int_t i = GetIndexFromLabel(lab);
158  return GetAcceptCluster(i);
159 }
160 
166 {
167  const Int_t n = GetNEntries();
168  AliVCluster *c = 0;
169  do {
170  fCurrentID++;
171  if (fCurrentID >= n) break;
173  } while (!c);
174 
175  return c;
176 }
177 
183 {
184  const Int_t n = GetNEntries();
185  AliVCluster *c = 0;
186  do {
187  fCurrentID++;
188  if (fCurrentID >= n) break;
189  c = GetCluster(fCurrentID);
190  } while (!c);
191 
192  return c;
193 }
194 
195 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
196 {
197  if (mass < 0) mass = 0;
198 
199  Double_t energy = 0;
200 
201  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
202  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
203  }
204  else {
205  energy = vc->E();
206  }
207 
208  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
209 
210  Float_t pos[3];
211  vc->GetPosition(pos);
212 
213  pos[0]-=fVertex[0];
214  pos[1]-=fVertex[1];
215  pos[2]-=fVertex[2];
216 
217  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
218 
219  if (r > 1e-12) {
220  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
221  }
222  else {
223  AliInfo("Null cluster radius, momentum calculation not possible");
224  return kFALSE;
225  }
226 
227  return kTRUE;
228 }
229 
230 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
231 {
232  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
233 
234  if (vc) {
235  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
236  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
237  }
238  else {
239  vc->GetMomentum(mom, fVertex);
240  }
241  }
242  else {
243  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
244  return kFALSE;
245  }
246  return kTRUE;
247 }
248 
255 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
256 {
257  AliVCluster *vc = GetCluster(i);
258  return GetMomentum(mom, vc);
259 }
260 
266 Bool_t AliClusterContainer::GetNextMomentum(TLorentzVector &mom)
267 {
268  AliVCluster *vc = GetNextCluster();
269  return GetMomentum(mom, vc);
270 }
271 
278 Bool_t AliClusterContainer::GetAcceptMomentum(TLorentzVector &mom, Int_t i) const
279 {
280  AliVCluster *vc = GetAcceptCluster(i);
281  return GetMomentum(mom, vc);
282 }
283 
290 {
291  AliVCluster *vc = GetNextAcceptCluster();
292  return GetMomentum(mom, vc);
293 }
294 
295 Bool_t AliClusterContainer::AcceptCluster(Int_t i, UInt_t &rejectionReason) const
296 {
297  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
298  if (!r) return kFALSE;
299 
300  AliTLorentzVector mom;
301  GetMomentum(mom, i);
302 
303  return ApplyKinematicCuts(mom, rejectionReason);
304 }
305 
306 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
307 {
308  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
309  if (!r) return kFALSE;
310 
311  AliTLorentzVector mom;
312  GetMomentum(mom, clus);
313 
314  return ApplyKinematicCuts(mom, rejectionReason);
315 }
316 
322 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
323 {
324  if (!clus) {
325  rejectionReason |= kNullObject;
326  return kFALSE;
327  }
328 
329  if (!clus->IsEMCAL()) {
330  rejectionReason |= kIsEMCalCut;
331  return kFALSE;
332  }
333 
334  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
335  rejectionReason |= kBitMapCut;
336  return kFALSE;
337  }
338 
339  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMinMCLabel) {
340  rejectionReason |= kMCLabelCut;
341  return kFALSE;
342  }
343 
344  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMaxMCLabel) {
345  rejectionReason |= kMCLabelCut;
346  return kFALSE;
347  }
348 
349  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
350  rejectionReason |= kTimeCut;
351  return kFALSE;
352  }
353 
354  if (fExoticCut && clus->GetIsExotic()) {
355  rejectionReason |= kExoticCut;
356  return kFALSE;
357  }
358 
359  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
360  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
361  rejectionReason |= kEnergyCut;
362  return kFALSE;
363  }
364  }
365 
366  return kTRUE;
367 }
368 
374 {
375  UInt_t rejectionReason = 0;
376  Int_t nClus = 0;
377  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
378  AliVCluster *clust = this->GetCluster(iclust);
379  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
380  }
381  return nClus;
382 }
383 
385 {
386  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
387  return fUserDefEnergyCut[t];
388  }
389  else {
390  return fMinE;
391  }
392 }
393 
395 {
396  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
397  fUserDefEnergyCut[t] = cut;
398  }
399  else {
400  fMinE = cut;
401  }
402 }
403 
404 const char* AliClusterContainer::GetTitle() const
405 {
406  static TString clusterString;
407 
409 
410  if (Ecut == 0) {
411  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
412  }
413  else if (Ecut < 1.0) {
414  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
415  }
416  else {
417  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
418  }
419 
420  return clusterString.Data();
421 }
422 
429  return all_iterator(this, 0, true);
430 }
431 
438  return all_iterator(this, this->GetNClusters(), true);
439 }
440 
447  return all_iterator(this, this->GetNClusters()-1, false);
448 }
449 
456  return all_iterator(this, -1, false);
457 }
458 
465  return accept_iterator(this, 0, true);
466 }
467 
474  return accept_iterator(this, this->GetNAcceptedClusters(), true);
475 }
476 
483  return accept_iterator(this, this->GetNAcceptedClusters()-1, false);
484 }
485 
492  return accept_iterator(this, -1, false);
493 }
494 
495 
501 
518  const AliClusterContainer *cont,
519  int start,
520  bool forward):
521  fkContainer(cont),
522  fAcceptIndices(),
523  fCurrentPos(start),
524  fForward(forward)
525 {
527  UInt_t rejection = 0;
528  Int_t mappos = 0;
529  for(int i = 0; i < fkContainer->GetNEntries(); i++){
530  if(fkContainer->AcceptCluster(i, rejection)){
531  fAcceptIndices[mappos++] = i;
532  }
533  }
534 }
535 
542  fkContainer(other.fkContainer),
543  fAcceptIndices(other.fAcceptIndices),
544  fCurrentPos(other.fCurrentPos),
545  fForward(other.fForward)
546 {
547 }
548 
556  if(this != &other){
557  fkContainer = other.fkContainer;
558  fAcceptIndices = other.fAcceptIndices;
559  fCurrentPos = other.fCurrentPos;
560  fForward = other.fForward;
561  }
562  return *this;
563 }
564 
572  return fCurrentPos != other.fCurrentPos;
573 }
574 
581  if(fForward){
582  fCurrentPos++;
583  } else {
584  fCurrentPos--;
585  }
586  return *this;
587 }
588 
597  this->operator++();
598  return result;
599 }
600 
607  if(fForward){
608  fCurrentPos--;
609  } else {
610  fCurrentPos++;
611  }
612  return *this;
613 }
614 
623  this->operator++();
624  return result;
625 }
626 
633  if(fCurrentPos < 0 || fCurrentPos >= fAcceptIndices.GetSize()){
634  return NULL;
635  }
636  return fkContainer->GetCluster(fAcceptIndices[fCurrentPos]);
637 }
638 
644 
659  const AliClusterContainer *cont,
660  int startpos,
661  bool forward):
662  fkContainer(cont),
663  fCurrentPos(startpos),
664  fForward(forward)
665 {
666 }
667 
674  fkContainer(other.fkContainer),
675  fCurrentPos(other.fCurrentPos),
676  fForward(other.fForward)
677 {
678 }
679 
687  if(&other != this){
688  fkContainer = other.fkContainer;
689  fCurrentPos = other.fCurrentPos;
690  fForward = other.fForward;
691  }
692  return *this;
693 }
694 
702  return other.fCurrentPos != fCurrentPos;
703 }
704 
711  if(fForward){
712  fCurrentPos++;
713  } else {
714  fCurrentPos--;
715  }
716  return *this;
717 }
718 
726  all_iterator tmp(*this);
727  operator++();
728  return tmp;
729 }
730 
737  if(fForward){
738  fCurrentPos--;
739  } else {
740  fCurrentPos++;
741  }
742  return *this;
743 }
744 
752  all_iterator tmp(*this);
753  operator--();
754  return tmp;
755 }
756 
763  if(fCurrentPos >= 0 && fCurrentPos < fkContainer->GetNClusters())
764  return fkContainer->GetCluster(fCurrentPos);
765  return NULL;
766 }
767 
768 /******************************************
769  * Unit tests *
770  ******************************************/
771 
772 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
773  std::vector<AliVCluster *> reference, variation;
774  AliVCluster *test = NULL;
775  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
776  test = cont->GetCluster(iclust);
777  if(!iteratorType){
778  UInt_t rejectionReason = 0;
779  if(!cont->AcceptCluster(test, rejectionReason)) continue;
780  }
781  reference.push_back(test);
782  }
783 
784  if(!iteratorType){
785  // test accept iterator
786  for(AliClusterContainer::accept_iterator iter = cont->accept_begin(); iter != cont->accept_end(); ++iter){
787  variation.push_back(*iter);
788  }
789  } else {
790  // test all iterator
791  for(AliClusterContainer::all_iterator iter = cont->begin(); iter != cont->end(); ++iter){
792  variation.push_back(*iter);
793  }
794  }
795 
796  if(verbose){
797  // Printing cluster addresses:
798  if(reference.size() < 30){
799  std::cout << "Clusters in reference container: " << std::endl;
800  std::cout << "===========================================" << std::endl;
801  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
802  std::cout << "Address: " << *refit << std::endl;
803  }
804  }
805  if(variation.size() < 30){
806  std::cout << "Clusters in test container: " << std::endl;
807  std::cout << "===========================================" << std::endl;
808  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
809  std::cout << "Address: " << *varit << std::endl;
810  }
811  }
812  }
813 
814  int testresult = 0;
815  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
816  bool failure = false;
817  // first test: cleck if all clusters are found by the iterator
818  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
819  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
820  if(verbose)
821  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
822  failure = true;
823  break;
824  }
825  }
826  if(!failure){
827  // second test: check if there are no extra clusters found
828  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
829  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
830  if(verbose)
831  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
832  failure = true;
833  break;
834  }
835  }
836  if(failure) testresult = 2;
837  } else {
838  testresult = 1;
839  }
840 
841  if(verbose){
842  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
843  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
844  std::cout << "Test result: " << testresult << std::endl;
845  std::cout << "-----------------------------------------------------------------------" << std::endl;
846  }
847  return testresult;
848 }
AliVCluster * GetAcceptClusterWithLabel(Int_t lab) const
bool fForward
Direction, expressed in forward direction.
bool fForward
Direction, expressed in forward direction.
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
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
stl iterator over accepted clusters
const char * GetTitle() const
all_iterator begin() const
accept_iterator accept_rbegin() const
Double_t fMinE
Max. cut on particle energy.
Bool_t GetNextMomentum(TLorentzVector &mom)
Int_t GetIndexFromLabel(Int_t lab) const
Double_t mass
bool operator!=(const accept_iterator &other) const
accept_iterator accept_end() const
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) ...
const AliClusterContainer * fkContainer
Container iterated over.
accept_iterator accept_begin() const
all_iterator end() const
int fCurrentPos
Current position inside the container.
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
accept_iterator & operator=(const accept_iterator &other)
enum AliVCluster::VCluUserDefEnergy_t VCluUserDefEnergy_t
Cluster is exotic cluster.
Energy below threshold.
stl iterator over all clusters
all_iterator & operator=(const all_iterator &other)
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
int fCurrentPos
Current position inside the container.
all_iterator rend() const
Double_t fClusTimeCutUp
up time cut for clusters
AliVCluster * GetCluster(Int_t i) const
void SetClassName(const char *clname)
Bool_t fExoticCut
reject clusters marked as "exotic"
AliVCluster * GetNextCluster()
bool operator!=(const all_iterator &other) const
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)
energy
Int_t fMinMCLabel
minimum MC label
accept_iterator accept_rend() const
const TString & GetArrayName() const
Double_t fMassHypothesis
if < 0 it will use a PID mass when available
Cluster not in the EMCAL.
Double_t fVertex[3]
! event vertex array
Bool_t GetAcceptMomentum(TLorentzVector &mom, Int_t i) const
AliVCluster * GetClusterWithLabel(Int_t lab) const
int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose)
TArrayI fAcceptIndices
Indices of accepted clusters.
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
TClonesArray * fClArray
! Pointer to array in input event
Int_t GetNEntries() const
Int_t fDefaultClusterEnergy
default cluster energy: -1 for clus->E(); otherwise clus->GetUserDefEnergy(fDefaultClusterEnergy) ...
const AliClusterContainer * fkContainer
Container iterated over.
all_iterator rbegin() const
Int_t fCurrentID
! current ID for automatic loops
Container structure for EMCAL clusters.
AliVCluster * GetNextAcceptCluster()
void ResetCurrentID(Int_t i=-1)