AliPhysics  a9863a5 (a9863a5)
 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 
423 
424 /******************************************
425  * Unit tests *
426  ******************************************/
427 
428 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
429  std::vector<AliVCluster *> reference, variation;
430  AliVCluster *test = NULL;
431  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
432  test = cont->GetCluster(iclust);
433  if(!iteratorType){
434  UInt_t rejectionReason = 0;
435  if(!cont->AcceptCluster(test, rejectionReason)) continue;
436  }
437  reference.push_back(test);
438  }
439 
440  if(!iteratorType){
441  // test accept iterator
442  for(AliClusterIterableContainer::iterator iter = cont->accept_begin(); iter != cont->accept_end(); ++iter){
443  variation.push_back(*iter);
444  }
445  } else {
446  // test all iterator
447  for(AliClusterIterableContainer::iterator iter = cont->begin(); iter != cont->end(); ++iter){
448  variation.push_back(*iter);
449  }
450  }
451 
452  if(verbose){
453  // Printing cluster addresses:
454  if(reference.size() < 30){
455  std::cout << "Clusters in reference container: " << std::endl;
456  std::cout << "===========================================" << std::endl;
457  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
458  std::cout << "Address: " << *refit << std::endl;
459  }
460  }
461  if(variation.size() < 30){
462  std::cout << "Clusters in test container: " << std::endl;
463  std::cout << "===========================================" << std::endl;
464  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
465  std::cout << "Address: " << *varit << std::endl;
466  }
467  }
468  }
469 
470  int testresult = 0;
471  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
472  bool failure = false;
473  // first test: cleck if all clusters are found by the iterator
474  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
475  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
476  if(verbose)
477  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
478  failure = true;
479  break;
480  }
481  }
482  if(!failure){
483  // second test: check if there are no extra clusters found
484  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
485  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
486  if(verbose)
487  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
488  failure = true;
489  break;
490  }
491  }
492  if(failure) testresult = 2;
493  } else {
494  testresult = 1;
495  }
496 
497  if(verbose){
498  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
499  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
500  std::cout << "Test result: " << testresult << std::endl;
501  std::cout << "-----------------------------------------------------------------------" << std::endl;
502  }
503  return testresult;
504 }
AliVCluster * GetAcceptClusterWithLabel(Int_t lab) const
AliClusterIterableContainer::iterator begin() const
virtual Bool_t ApplyClusterCuts(const AliVCluster *clus, UInt_t &rejectionReason) const
Bool_t GetNextAcceptMomentum(TLorentzVector &mom)
AliClusterIterableContainer::iterator end() const
virtual Bool_t ApplyKinematicCuts(const AliTLorentzVector &mom, UInt_t &rejectionReason) const
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
const char * GetTitle() const
Double_t fMinE
Max. cut on particle energy.
AliClusterIterableContainer::iterator accept_end() const
Bool_t GetNextMomentum(TLorentzVector &mom)
Int_t GetIndexFromLabel(Int_t lab) const
Double_t mass
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) ...
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
bidirectional stl iterator over the EMCAL iterable container
Int_t GetDefaultClusterEnergy() const
enum AliVCluster::VCluUserDefEnergy_t VCluUserDefEnergy_t
Cluster is exotic cluster.
Energy below threshold.
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
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_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
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
AliClusterIterableContainer::iterator accept_begin() const
AliVCluster * GetClusterWithLabel(Int_t lab) const
int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose)
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) ...
Int_t fCurrentID
! current ID for automatic loops
Container structure for EMCAL clusters.
AliVCluster * GetNextAcceptCluster()
void ResetCurrentID(Int_t i=-1)