AliPhysics  3337bb0 (3337bb0)
 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 "AliAODEvent.h"
21 #include "AliESDEvent.h"
22 #include "AliVEvent.h"
23 #include "AliLog.h"
24 #include "AliTLorentzVector.h"
25 
26 #include "AliClusterContainer.h"
27 
29 ClassImp(AliClusterContainer);
31 
32 // Properly instantiate the object
34 
39  AliEmcalContainer(),
40  fClusTimeCutLow(-10),
41  fClusTimeCutUp(10),
42  fExoticCut(kTRUE),
43  fDefaultClusterEnergy(-1),
44  fIncludePHOS(kFALSE),
45  fIncludePHOSonly(kFALSE),
46  fPhosMinNcells(0),
47  fPhosMinM02(0)
48 {
49  fBaseClassName = "AliVCluster";
50  SetClassName("AliVCluster");
51 
52  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
53  fUserDefEnergyCut[i] = 0.;
54  }
55 }
56 
62  AliEmcalContainer(name),
63  fClusTimeCutLow(-10),
64  fClusTimeCutUp(10),
65  fExoticCut(kTRUE),
66  fDefaultClusterEnergy(-1),
67  fIncludePHOS(kFALSE),
68  fIncludePHOSonly(kFALSE),
69  fPhosMinNcells(0),
70  fPhosMinM02(0)
71 {
72  fBaseClassName = "AliVCluster";
73  SetClassName("AliVCluster");
74 
75  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
76  fUserDefEnergyCut[i] = 0.;
77  }
78 }
79 
85 AliVCluster* AliClusterContainer::GetLeadingCluster(const char* opt)
86 {
87  TString option(opt);
88  option.ToLower();
89 
90  double (AliTLorentzVector::*momentum)() const = 0;
91 
92  if (option.Contains("e")) {
93  momentum = &AliTLorentzVector::E;
94  }
95  else {
96  momentum = &AliTLorentzVector::Et;
97  }
98 
100 
101  for (auto cluster : accepted_momentum()) {
102  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
103  clusterMax = cluster;
104  }
105  }
106 
107  return clusterMax.second;
108 }
109 
116 {
117  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
118  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
119  return vp;
120 
121 }
122 
129 {
130  AliVCluster *vc = GetCluster(i);
131  if (!vc) return 0;
132 
133  UInt_t rejectionReason = 0;
134  if (AcceptCluster(vc, rejectionReason))
135  return vc;
136  else {
137  AliDebug(2,"Cluster not accepted.");
138  return 0;
139  }
140 }
141 
148 {
149  Int_t i = GetIndexFromLabel(lab);
150  return GetCluster(i);
151 }
152 
159 {
160  Int_t i = GetIndexFromLabel(lab);
161  return GetAcceptCluster(i);
162 }
163 
169 {
170  const Int_t n = GetNEntries();
171  AliVCluster *c = 0;
172  do {
173  fCurrentID++;
174  if (fCurrentID >= n) break;
175  c = GetAcceptCluster(fCurrentID);
176  } while (!c);
177 
178  return c;
179 }
180 
186 {
187  const Int_t n = GetNEntries();
188  AliVCluster *c = 0;
189  do {
190  fCurrentID++;
191  if (fCurrentID >= n) break;
192  c = GetCluster(fCurrentID);
193  } while (!c);
194 
195  return c;
196 }
197 
198 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
199 {
200  if (mass < 0) mass = 0;
201 
202  Double_t energy = 0;
203 
204  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
205  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
206  }
207  else {
208  energy = vc->E();
209  }
210 
211  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
212 
213  Float_t pos[3];
214  vc->GetPosition(pos);
215 
216  pos[0]-=fVertex[0];
217  pos[1]-=fVertex[1];
218  pos[2]-=fVertex[2];
219 
220  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
221 
222  if (r > 1e-12) {
223  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
224  }
225  else {
226  AliInfo("Null cluster radius, momentum calculation not possible");
227  return kFALSE;
228  }
229 
230  return kTRUE;
231 }
232 
233 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
234 {
235  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
236 
237  if (vc) {
238  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
239  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
240  }
241  else {
242  vc->GetMomentum(mom, fVertex);
243  }
244  }
245  else {
246  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
247  return kFALSE;
248  }
249  return kTRUE;
250 }
251 
258 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
259 {
260  AliVCluster *vc = GetCluster(i);
261  return GetMomentum(mom, vc);
262 }
263 
270 {
271  AliVCluster *vc = GetNextCluster();
272  return GetMomentum(mom, vc);
273 }
274 
282 {
283  AliVCluster *vc = GetAcceptCluster(i);
284  return GetMomentum(mom, vc);
285 }
286 
293 {
294  AliVCluster *vc = GetNextAcceptCluster();
295  return GetMomentum(mom, vc);
296 }
297 
299 {
300  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
301  if (!r) return kFALSE;
302 
303  AliTLorentzVector mom;
304  GetMomentum(mom, i);
305 
306  return ApplyKinematicCuts(mom, rejectionReason);
307 }
308 
309 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
310 {
311  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
312  if (!r) return kFALSE;
313 
314  AliTLorentzVector mom;
315  GetMomentum(mom, clus);
316 
317  return ApplyKinematicCuts(mom, rejectionReason);
318 }
319 
325 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
326 {
327  if (!clus) {
328  rejectionReason |= kNullObject;
329  return kFALSE;
330  }
331 
332  Bool_t bInAcceptance = clus->IsEMCAL();
333  if (fIncludePHOS) {
334  bInAcceptance = clus->IsEMCAL() || (clus->GetType() == AliVCluster::kPHOSNeutral);
335  }
336  if (fIncludePHOSonly) {
337  bInAcceptance = (clus->GetType() == AliVCluster::kPHOSNeutral);
338  }
339  if (!bInAcceptance) {
340  rejectionReason |= kIsEMCalCut;
341  return kFALSE;
342  }
343 
344  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
345  rejectionReason |= kBitMapCut;
346  return kFALSE;
347  }
348 
349  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
350  rejectionReason |= kMCLabelCut;
351  return kFALSE;
352  }
353 
354  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
355  rejectionReason |= kMCLabelCut;
356  return kFALSE;
357  }
358 
359  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
360  rejectionReason |= kTimeCut;
361  return kFALSE;
362  }
363 
364  if (fExoticCut && clus->GetIsExotic()) {
365  rejectionReason |= kExoticCut;
366  return kFALSE;
367  }
368 
369  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
370  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
371  rejectionReason |= kEnergyCut;
372  return kFALSE;
373  }
374  }
375 
377  if (clus->GetType() == AliVCluster::kPHOSNeutral) {
378  if (clus->GetNCells() < fPhosMinNcells) {
379  rejectionReason |= kExoticCut;
380  return kFALSE;
381  }
382 
383  if (clus->GetM02() < fPhosMinM02) {
384  rejectionReason |= kExoticCut;
385  return kFALSE;
386  }
387  }
388  }
389 
390  return kTRUE;
391 }
392 
398 {
399  UInt_t rejectionReason = 0;
400  Int_t nClus = 0;
401  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
402  AliVCluster *clust = this->GetCluster(iclust);
403  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
404  }
405  return nClus;
406 }
407 
414 {
415  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
416  return fUserDefEnergyCut[t];
417  }
418  else {
419  return fMinE;
420  }
421 }
422 
429 {
430  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
431  fUserDefEnergyCut[t] = cut;
432  }
433  else {
434  fMinE = cut;
435  }
436 }
437 
446 void AliClusterContainer::SetArray(const AliVEvent * event)
447 {
448  AliEmcalContainer::SetArray(event);
449 
450  // Register TClonesArray in index map
452 }
453 
460  return AliClusterIterableContainer(this, false);
461 }
462 
469  return AliClusterIterableContainer(this, true);
470 }
471 
478  return AliClusterIterableMomentumContainer(this, false);
479 }
480 
487  return AliClusterIterableMomentumContainer(this, true);
488 }
489 
490 const char* AliClusterContainer::GetTitle() const
491 {
492  static TString clusterString;
493 
495 
496  if (Ecut == 0) {
497  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
498  }
499  else if (Ecut < 1.0) {
500  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
501  }
502  else {
503  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
504  }
505 
506  return clusterString.Data();
507 }
508 
509 TString AliClusterContainer::GetDefaultArrayName(const AliVEvent * const ev) const {
510  if(ev->IsA() == AliAODEvent::Class()) return "caloClusters";
511  else if(ev->IsA() == AliESDEvent::Class()) return "CaloClusters";
512  else return "";
513 }
514 
515 
516 /******************************************
517  * Unit tests *
518  ******************************************/
519 
520 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
521  std::vector<AliVCluster *> reference, variation;
522  AliVCluster *test = NULL;
523  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
524  test = cont->GetCluster(iclust);
525  if(!iteratorType){
526  UInt_t rejectionReason = 0;
527  if(!cont->AcceptCluster(test, rejectionReason)) continue;
528  }
529  reference.push_back(test);
530  }
531 
532  if(!iteratorType){
533  // test accept iterator
534  for(auto cluster : cont->accepted()){
535  variation.push_back(cluster);
536  }
537  } else {
538  // test all iterator
539  for(auto cluster : cont->all()){
540  variation.push_back(cluster);
541  }
542  }
543 
544  if(verbose){
545  // Printing cluster addresses:
546  if(reference.size() < 30){
547  std::cout << "Clusters in reference container: " << std::endl;
548  std::cout << "===========================================" << std::endl;
549  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
550  std::cout << "Address: " << *refit << std::endl;
551  }
552  }
553  if(variation.size() < 30){
554  std::cout << "Clusters in test container: " << std::endl;
555  std::cout << "===========================================" << std::endl;
556  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
557  std::cout << "Address: " << *varit << std::endl;
558  }
559  }
560  }
561 
562  int testresult = 0;
563  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
564  bool failure = false;
565  // first test: cleck if all clusters are found by the iterator
566  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
567  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
568  if(verbose)
569  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
570  failure = true;
571  break;
572  }
573  }
574  if(!failure){
575  // second test: check if there are no extra clusters found
576  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
577  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
578  if(verbose)
579  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
580  failure = true;
581  break;
582  }
583  }
584  if(failure) testresult = 2;
585  } else {
586  testresult = 1;
587  }
588 
589  if(verbose){
590  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
591  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
592  std::cout << "Test result: " << testresult << std::endl;
593  std::cout << "-----------------------------------------------------------------------" << std::endl;
594  }
595  return testresult;
596 }
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
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