AliPhysics  6eb37c2 (6eb37c2)
 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 {
52  fBaseClassName = "AliVCluster";
53  SetClassName("AliVCluster");
54 
55  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
56  fUserDefEnergyCut[i] = 0.;
57  }
58 }
59 
65  AliEmcalContainer(name),
66  fClusTimeCutLow(-10),
67  fClusTimeCutUp(10),
68  fExoticCut(kTRUE),
69  fDefaultClusterEnergy(-1),
70  fIncludePHOS(kFALSE),
71  fIncludePHOSonly(kFALSE),
72  fPhosMinNcells(0),
73  fPhosMinM02(0),
74  fEmcalMinM02(DBL_MIN),
75  fEmcalMaxM02(DBL_MAX)
76 {
77  fBaseClassName = "AliVCluster";
78  SetClassName("AliVCluster");
79 
80  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
81  fUserDefEnergyCut[i] = 0.;
82  }
83 }
84 
90 AliVCluster* AliClusterContainer::GetLeadingCluster(const char* opt)
91 {
92  TString option(opt);
93  option.ToLower();
94 
95  double (AliTLorentzVector::*momentum)() const = 0;
96 
97  if (option.Contains("e")) {
98  momentum = &AliTLorentzVector::E;
99  }
100  else {
101  momentum = &AliTLorentzVector::Et;
102  }
103 
105 
106  for (auto cluster : accepted_momentum()) {
107  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
108  clusterMax = cluster;
109  }
110  }
111 
112  return clusterMax.second;
113 }
114 
121 {
122  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
123  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
124  return vp;
125 
126 }
127 
134 {
135  AliVCluster *vc = GetCluster(i);
136  if (!vc) return 0;
137 
138  UInt_t rejectionReason = 0;
139  if (AcceptCluster(vc, rejectionReason))
140  return vc;
141  else {
142  AliDebug(2,"Cluster not accepted.");
143  return 0;
144  }
145 }
146 
153 {
154  Int_t i = GetIndexFromLabel(lab);
155  return GetCluster(i);
156 }
157 
164 {
165  Int_t i = GetIndexFromLabel(lab);
166  return GetAcceptCluster(i);
167 }
168 
174 {
175  const Int_t n = GetNEntries();
176  AliVCluster *c = 0;
177  do {
178  fCurrentID++;
179  if (fCurrentID >= n) break;
180  c = GetAcceptCluster(fCurrentID);
181  } while (!c);
182 
183  return c;
184 }
185 
191 {
192  const Int_t n = GetNEntries();
193  AliVCluster *c = 0;
194  do {
195  fCurrentID++;
196  if (fCurrentID >= n) break;
197  c = GetCluster(fCurrentID);
198  } while (!c);
199 
200  return c;
201 }
202 
203 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
204 {
205  if (mass < 0) mass = 0;
206 
207  Double_t energy = 0;
208 
209  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
210  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
211  }
212  else {
213  energy = vc->E();
214  }
215 
216  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
217 
218  Float_t pos[3];
219  vc->GetPosition(pos);
220 
221  pos[0]-=fVertex[0];
222  pos[1]-=fVertex[1];
223  pos[2]-=fVertex[2];
224 
225  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
226 
227  if (r > 1e-12) {
228  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
229  }
230  else {
231  AliInfo("Null cluster radius, momentum calculation not possible");
232  return kFALSE;
233  }
234 
235  return kTRUE;
236 }
237 
238 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
239 {
240  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
241 
242  if (vc) {
243  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
244  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
245  }
246  else {
247  vc->GetMomentum(mom, fVertex);
248  }
249  }
250  else {
251  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
252  return kFALSE;
253  }
254  return kTRUE;
255 }
256 
263 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
264 {
265  AliVCluster *vc = GetCluster(i);
266  return GetMomentum(mom, vc);
267 }
268 
275 {
276  AliVCluster *vc = GetNextCluster();
277  return GetMomentum(mom, vc);
278 }
279 
287 {
288  AliVCluster *vc = GetAcceptCluster(i);
289  return GetMomentum(mom, vc);
290 }
291 
298 {
299  AliVCluster *vc = GetNextAcceptCluster();
300  return GetMomentum(mom, vc);
301 }
302 
304 {
305  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
306  if (!r) return kFALSE;
307 
308  AliTLorentzVector mom;
309  GetMomentum(mom, i);
310 
311  return ApplyKinematicCuts(mom, rejectionReason);
312 }
313 
314 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
315 {
316  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
317  if (!r) return kFALSE;
318 
319  AliTLorentzVector mom;
320  GetMomentum(mom, clus);
321 
322  return ApplyKinematicCuts(mom, rejectionReason);
323 }
324 
330 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
331 {
332  if (!clus) {
333  rejectionReason |= kNullObject;
334  return kFALSE;
335  }
336 
337  Bool_t bInAcceptance = clus->IsEMCAL();
338  if (fIncludePHOS) {
339  bInAcceptance = clus->IsEMCAL() || (clus->GetType() == AliVCluster::kPHOSNeutral);
340  }
341  if (fIncludePHOSonly) {
342  bInAcceptance = (clus->GetType() == AliVCluster::kPHOSNeutral);
343  }
344  if (!bInAcceptance) {
345  rejectionReason |= kIsEMCalCut;
346  return kFALSE;
347  }
348 
349  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
350  rejectionReason |= kBitMapCut;
351  return kFALSE;
352  }
353 
354  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
355  rejectionReason |= kMCLabelCut;
356  return kFALSE;
357  }
358 
359  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
360  rejectionReason |= kMCLabelCut;
361  return kFALSE;
362  }
363 
364  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
365  rejectionReason |= kTimeCut;
366  return kFALSE;
367  }
368 
369  if (fExoticCut && clus->GetIsExotic()) {
370  rejectionReason |= kExoticCut;
371  return kFALSE;
372  }
373 
374  if(clus->GetM02() < fEmcalMinM02 || clus->GetM02() > fEmcalMaxM02) {
375  rejectionReason |= kExoticCut; // Not really true, but there is a lack of leftover bits
376  return kFALSE;
377  }
378 
379  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
380  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
381  rejectionReason |= kEnergyCut;
382  return kFALSE;
383  }
384  }
385 
387  if (clus->GetType() == AliVCluster::kPHOSNeutral) {
388  if (clus->GetNCells() < fPhosMinNcells) {
389  rejectionReason |= kExoticCut;
390  return kFALSE;
391  }
392 
393  if (clus->GetM02() < fPhosMinM02) {
394  rejectionReason |= kExoticCut;
395  return kFALSE;
396  }
397  }
398  }
399 
400  return kTRUE;
401 }
402 
408 {
409  UInt_t rejectionReason = 0;
410  Int_t nClus = 0;
411  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
412  AliVCluster *clust = this->GetCluster(iclust);
413  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
414  }
415  return nClus;
416 }
417 
424 {
425  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
426  return fUserDefEnergyCut[t];
427  }
428  else {
429  return fMinE;
430  }
431 }
432 
439 {
440  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
441  fUserDefEnergyCut[t] = cut;
442  }
443  else {
444  fMinE = cut;
445  }
446 }
447 
456 void AliClusterContainer::SetArray(const AliVEvent * event)
457 {
458  AliEmcalContainer::SetArray(event);
459 
460  // Register TClonesArray in index map
462 }
463 
470  return AliClusterIterableContainer(this, false);
471 }
472 
479  return AliClusterIterableContainer(this, true);
480 }
481 
488  return AliClusterIterableMomentumContainer(this, false);
489 }
490 
497  return AliClusterIterableMomentumContainer(this, true);
498 }
499 
500 const char* AliClusterContainer::GetTitle() const
501 {
502  static TString clusterString;
503 
505 
506  if (Ecut == 0) {
507  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
508  }
509  else if (Ecut < 1.0) {
510  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
511  }
512  else {
513  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
514  }
515 
516  return clusterString.Data();
517 }
518 
519 TString AliClusterContainer::GetDefaultArrayName(const AliVEvent * const ev) const {
520  if(ev->IsA() == AliAODEvent::Class()) return "caloClusters";
521  else if(ev->IsA() == AliESDEvent::Class()) return "CaloClusters";
522  else return "";
523 }
524 
525 
526 /******************************************
527  * Unit tests *
528  ******************************************/
529 
530 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
531  std::vector<AliVCluster *> reference, variation;
532  AliVCluster *test = NULL;
533  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
534  test = cont->GetCluster(iclust);
535  if(!iteratorType){
536  UInt_t rejectionReason = 0;
537  if(!cont->AcceptCluster(test, rejectionReason)) continue;
538  }
539  reference.push_back(test);
540  }
541 
542  if(!iteratorType){
543  // test accept iterator
544  for(auto cluster : cont->accepted()){
545  variation.push_back(cluster);
546  }
547  } else {
548  // test all iterator
549  for(auto cluster : cont->all()){
550  variation.push_back(cluster);
551  }
552  }
553 
554  if(verbose){
555  // Printing cluster addresses:
556  if(reference.size() < 30){
557  std::cout << "Clusters in reference container: " << std::endl;
558  std::cout << "===========================================" << std::endl;
559  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
560  std::cout << "Address: " << *refit << std::endl;
561  }
562  }
563  if(variation.size() < 30){
564  std::cout << "Clusters in test container: " << std::endl;
565  std::cout << "===========================================" << std::endl;
566  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
567  std::cout << "Address: " << *varit << std::endl;
568  }
569  }
570  }
571 
572  int testresult = 0;
573  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
574  bool failure = false;
575  // first test: cleck if all clusters are found by the iterator
576  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
577  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
578  if(verbose)
579  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
580  failure = true;
581  break;
582  }
583  }
584  if(!failure){
585  // second test: check if there are no extra clusters found
586  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
587  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
588  if(verbose)
589  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
590  failure = true;
591  break;
592  }
593  }
594  if(failure) testresult = 2;
595  } else {
596  testresult = 1;
597  }
598 
599  if(verbose){
600  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
601  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
602  std::cout << "Test result: " << testresult << std::endl;
603  std::cout << "-----------------------------------------------------------------------" << std::endl;
604  }
605  return testresult;
606 }
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
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