AliPhysics  db95e02 (db95e02)
 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 
34  AliEmcalContainer(),
35  fClusTimeCutLow(-10),
36  fClusTimeCutUp(10),
37  fExoticCut(kTRUE),
38  fDefaultClusterEnergy(-1),
39  fIncludePHOS(kFALSE)
40 {
41  fBaseClassName = "AliVCluster";
42  SetClassName("AliVCluster");
43 
44  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
45  fUserDefEnergyCut[i] = 0.;
46  }
47 }
48 
54  AliEmcalContainer(name),
55  fClusTimeCutLow(-10),
56  fClusTimeCutUp(10),
57  fExoticCut(kTRUE),
58  fDefaultClusterEnergy(-1),
59  fIncludePHOS(kFALSE)
60 {
61  fBaseClassName = "AliVCluster";
62  SetClassName("AliVCluster");
63 
64  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
65  fUserDefEnergyCut[i] = 0.;
66  }
67 }
68 
74 AliVCluster* AliClusterContainer::GetLeadingCluster(const char* opt)
75 {
76  TString option(opt);
77  option.ToLower();
78 
79  double (AliTLorentzVector::*momentum)() const = 0;
80 
81  if (option.Contains("e")) {
82  momentum = &AliTLorentzVector::E;
83  }
84  else {
85  momentum = &AliTLorentzVector::Et;
86  }
87 
89 
90  for (auto cluster : accepted_momentum()) {
91  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
92  clusterMax = cluster;
93  }
94  }
95 
96  return clusterMax.second;
97 }
98 
105 {
106  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
107  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
108  return vp;
109 
110 }
111 
118 {
119  AliVCluster *vc = GetCluster(i);
120  if (!vc) return 0;
121 
122  UInt_t rejectionReason = 0;
123  if (AcceptCluster(vc, rejectionReason))
124  return vc;
125  else {
126  AliDebug(2,"Cluster not accepted.");
127  return 0;
128  }
129 }
130 
137 {
138  Int_t i = GetIndexFromLabel(lab);
139  return GetCluster(i);
140 }
141 
148 {
149  Int_t i = GetIndexFromLabel(lab);
150  return GetAcceptCluster(i);
151 }
152 
158 {
159  const Int_t n = GetNEntries();
160  AliVCluster *c = 0;
161  do {
162  fCurrentID++;
163  if (fCurrentID >= n) break;
164  c = GetAcceptCluster(fCurrentID);
165  } while (!c);
166 
167  return c;
168 }
169 
175 {
176  const Int_t n = GetNEntries();
177  AliVCluster *c = 0;
178  do {
179  fCurrentID++;
180  if (fCurrentID >= n) break;
181  c = GetCluster(fCurrentID);
182  } while (!c);
183 
184  return c;
185 }
186 
187 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
188 {
189  if (mass < 0) mass = 0;
190 
191  Double_t energy = 0;
192 
193  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
194  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
195  }
196  else {
197  energy = vc->E();
198  }
199 
200  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
201 
202  Float_t pos[3];
203  vc->GetPosition(pos);
204 
205  pos[0]-=fVertex[0];
206  pos[1]-=fVertex[1];
207  pos[2]-=fVertex[2];
208 
209  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
210 
211  if (r > 1e-12) {
212  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
213  }
214  else {
215  AliInfo("Null cluster radius, momentum calculation not possible");
216  return kFALSE;
217  }
218 
219  return kTRUE;
220 }
221 
222 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
223 {
224  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
225 
226  if (vc) {
227  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
228  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
229  }
230  else {
231  vc->GetMomentum(mom, fVertex);
232  }
233  }
234  else {
235  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
236  return kFALSE;
237  }
238  return kTRUE;
239 }
240 
247 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
248 {
249  AliVCluster *vc = GetCluster(i);
250  return GetMomentum(mom, vc);
251 }
252 
259 {
260  AliVCluster *vc = GetNextCluster();
261  return GetMomentum(mom, vc);
262 }
263 
271 {
272  AliVCluster *vc = GetAcceptCluster(i);
273  return GetMomentum(mom, vc);
274 }
275 
282 {
283  AliVCluster *vc = GetNextAcceptCluster();
284  return GetMomentum(mom, vc);
285 }
286 
288 {
289  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
290  if (!r) return kFALSE;
291 
292  AliTLorentzVector mom;
293  GetMomentum(mom, i);
294 
295  return ApplyKinematicCuts(mom, rejectionReason);
296 }
297 
298 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
299 {
300  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
301  if (!r) return kFALSE;
302 
303  AliTLorentzVector mom;
304  GetMomentum(mom, clus);
305 
306  return ApplyKinematicCuts(mom, rejectionReason);
307 }
308 
314 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
315 {
316  if (!clus) {
317  rejectionReason |= kNullObject;
318  return kFALSE;
319  }
320 
321  Bool_t bInAcceptance = clus->IsEMCAL();
322  if (fIncludePHOS) bInAcceptance = clus->IsEMCAL() || clus->IsPHOS();
323  if (!bInAcceptance) {
324  rejectionReason |= kIsEMCalCut;
325  return kFALSE;
326  }
327 
328  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
329  rejectionReason |= kBitMapCut;
330  return kFALSE;
331  }
332 
333  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
334  rejectionReason |= kMCLabelCut;
335  return kFALSE;
336  }
337 
338  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
339  rejectionReason |= kMCLabelCut;
340  return kFALSE;
341  }
342 
343  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
344  rejectionReason |= kTimeCut;
345  return kFALSE;
346  }
347 
348  if (fExoticCut && clus->GetIsExotic()) {
349  rejectionReason |= kExoticCut;
350  return kFALSE;
351  }
352 
353  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
354  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
355  rejectionReason |= kEnergyCut;
356  return kFALSE;
357  }
358  }
359 
360  return kTRUE;
361 }
362 
368 {
369  UInt_t rejectionReason = 0;
370  Int_t nClus = 0;
371  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
372  AliVCluster *clust = this->GetCluster(iclust);
373  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
374  }
375  return nClus;
376 }
377 
384 {
385  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
386  return fUserDefEnergyCut[t];
387  }
388  else {
389  return fMinE;
390  }
391 }
392 
399 {
400  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
401  fUserDefEnergyCut[t] = cut;
402  }
403  else {
404  fMinE = cut;
405  }
406 }
407 
414  return AliClusterIterableContainer(this, false);
415 }
416 
423  return AliClusterIterableContainer(this, true);
424 }
425 
432  return AliClusterIterableMomentumContainer(this, false);
433 }
434 
441  return AliClusterIterableMomentumContainer(this, true);
442 }
443 
444 const char* AliClusterContainer::GetTitle() const
445 {
446  static TString clusterString;
447 
449 
450  if (Ecut == 0) {
451  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
452  }
453  else if (Ecut < 1.0) {
454  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
455  }
456  else {
457  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
458  }
459 
460  return clusterString.Data();
461 }
462 
463 
464 /******************************************
465  * Unit tests *
466  ******************************************/
467 
468 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
469  std::vector<AliVCluster *> reference, variation;
470  AliVCluster *test = NULL;
471  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
472  test = cont->GetCluster(iclust);
473  if(!iteratorType){
474  UInt_t rejectionReason = 0;
475  if(!cont->AcceptCluster(test, rejectionReason)) continue;
476  }
477  reference.push_back(test);
478  }
479 
480  if(!iteratorType){
481  // test accept iterator
482  for(auto cluster : cont->accepted()){
483  variation.push_back(cluster);
484  }
485  } else {
486  // test all iterator
487  for(auto cluster : cont->all()){
488  variation.push_back(cluster);
489  }
490  }
491 
492  if(verbose){
493  // Printing cluster addresses:
494  if(reference.size() < 30){
495  std::cout << "Clusters in reference container: " << std::endl;
496  std::cout << "===========================================" << std::endl;
497  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
498  std::cout << "Address: " << *refit << std::endl;
499  }
500  }
501  if(variation.size() < 30){
502  std::cout << "Clusters in test container: " << std::endl;
503  std::cout << "===========================================" << std::endl;
504  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
505  std::cout << "Address: " << *varit << std::endl;
506  }
507  }
508  }
509 
510  int testresult = 0;
511  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
512  bool failure = false;
513  // first test: cleck if all clusters are found by the iterator
514  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
515  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
516  if(verbose)
517  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
518  failure = true;
519  break;
520  }
521  }
522  if(!failure){
523  // second test: check if there are no extra clusters found
524  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
525  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
526  if(verbose)
527  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
528  failure = true;
529  break;
530  }
531  }
532  if(failure) testresult = 2;
533  } else {
534  testresult = 1;
535  }
536 
537  if(verbose){
538  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
539  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
540  std::cout << "Test result: " << testresult << std::endl;
541  std::cout << "-----------------------------------------------------------------------" << std::endl;
542  }
543  return testresult;
544 }
AliVCluster * GetAcceptClusterWithLabel(Int_t lab) 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
Bool_t fIncludePHOS
whether or not to include 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)
energy
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
int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
Bool_t GetMomentum(TLorentzVector &mom, const AliVCluster *vc, Double_t mass) const
bool Bool_t
Definition: External.C:53
Int_t fDefaultClusterEnergy
default cluster energy: -1 for clus->E(); otherwise clus->GetUserDefEnergy(fDefaultClusterEnergy) ...
Container structure for EMCAL clusters.
AliVCluster * GetNextAcceptCluster()