AliPhysics  cdeda5a (cdeda5a)
 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 
30 // Properly instantiate the object
32 
37  AliEmcalContainer(),
38  fClusTimeCutLow(-10),
39  fClusTimeCutUp(10),
40  fExoticCut(kTRUE),
41  fDefaultClusterEnergy(-1),
42  fIncludePHOS(kFALSE)
43 {
44  fBaseClassName = "AliVCluster";
45  SetClassName("AliVCluster");
46 
47  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
48  fUserDefEnergyCut[i] = 0.;
49  }
50 }
51 
57  AliEmcalContainer(name),
58  fClusTimeCutLow(-10),
59  fClusTimeCutUp(10),
60  fExoticCut(kTRUE),
61  fDefaultClusterEnergy(-1),
62  fIncludePHOS(kFALSE)
63 {
64  fBaseClassName = "AliVCluster";
65  SetClassName("AliVCluster");
66 
67  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
68  fUserDefEnergyCut[i] = 0.;
69  }
70 }
71 
77 AliVCluster* AliClusterContainer::GetLeadingCluster(const char* opt)
78 {
79  TString option(opt);
80  option.ToLower();
81 
82  double (AliTLorentzVector::*momentum)() const = 0;
83 
84  if (option.Contains("e")) {
85  momentum = &AliTLorentzVector::E;
86  }
87  else {
88  momentum = &AliTLorentzVector::Et;
89  }
90 
92 
93  for (auto cluster : accepted_momentum()) {
94  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
95  clusterMax = cluster;
96  }
97  }
98 
99  return clusterMax.second;
100 }
101 
108 {
109  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
110  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
111  return vp;
112 
113 }
114 
121 {
122  AliVCluster *vc = GetCluster(i);
123  if (!vc) return 0;
124 
125  UInt_t rejectionReason = 0;
126  if (AcceptCluster(vc, rejectionReason))
127  return vc;
128  else {
129  AliDebug(2,"Cluster not accepted.");
130  return 0;
131  }
132 }
133 
140 {
141  Int_t i = GetIndexFromLabel(lab);
142  return GetCluster(i);
143 }
144 
151 {
152  Int_t i = GetIndexFromLabel(lab);
153  return GetAcceptCluster(i);
154 }
155 
161 {
162  const Int_t n = GetNEntries();
163  AliVCluster *c = 0;
164  do {
165  fCurrentID++;
166  if (fCurrentID >= n) break;
167  c = GetAcceptCluster(fCurrentID);
168  } while (!c);
169 
170  return c;
171 }
172 
178 {
179  const Int_t n = GetNEntries();
180  AliVCluster *c = 0;
181  do {
182  fCurrentID++;
183  if (fCurrentID >= n) break;
184  c = GetCluster(fCurrentID);
185  } while (!c);
186 
187  return c;
188 }
189 
190 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
191 {
192  if (mass < 0) mass = 0;
193 
194  Double_t energy = 0;
195 
196  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
197  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
198  }
199  else {
200  energy = vc->E();
201  }
202 
203  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
204 
205  Float_t pos[3];
206  vc->GetPosition(pos);
207 
208  pos[0]-=fVertex[0];
209  pos[1]-=fVertex[1];
210  pos[2]-=fVertex[2];
211 
212  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
213 
214  if (r > 1e-12) {
215  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
216  }
217  else {
218  AliInfo("Null cluster radius, momentum calculation not possible");
219  return kFALSE;
220  }
221 
222  return kTRUE;
223 }
224 
225 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
226 {
227  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
228 
229  if (vc) {
230  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
231  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
232  }
233  else {
234  vc->GetMomentum(mom, fVertex);
235  }
236  }
237  else {
238  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
239  return kFALSE;
240  }
241  return kTRUE;
242 }
243 
250 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
251 {
252  AliVCluster *vc = GetCluster(i);
253  return GetMomentum(mom, vc);
254 }
255 
262 {
263  AliVCluster *vc = GetNextCluster();
264  return GetMomentum(mom, vc);
265 }
266 
274 {
275  AliVCluster *vc = GetAcceptCluster(i);
276  return GetMomentum(mom, vc);
277 }
278 
285 {
286  AliVCluster *vc = GetNextAcceptCluster();
287  return GetMomentum(mom, vc);
288 }
289 
291 {
292  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
293  if (!r) return kFALSE;
294 
295  AliTLorentzVector mom;
296  GetMomentum(mom, i);
297 
298  return ApplyKinematicCuts(mom, rejectionReason);
299 }
300 
301 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
302 {
303  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
304  if (!r) return kFALSE;
305 
306  AliTLorentzVector mom;
307  GetMomentum(mom, clus);
308 
309  return ApplyKinematicCuts(mom, rejectionReason);
310 }
311 
317 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
318 {
319  if (!clus) {
320  rejectionReason |= kNullObject;
321  return kFALSE;
322  }
323 
324  Bool_t bInAcceptance = clus->IsEMCAL();
325  if (fIncludePHOS) bInAcceptance = clus->IsEMCAL() || clus->IsPHOS();
326  if (!bInAcceptance) {
327  rejectionReason |= kIsEMCalCut;
328  return kFALSE;
329  }
330 
331  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
332  rejectionReason |= kBitMapCut;
333  return kFALSE;
334  }
335 
336  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
337  rejectionReason |= kMCLabelCut;
338  return kFALSE;
339  }
340 
341  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
342  rejectionReason |= kMCLabelCut;
343  return kFALSE;
344  }
345 
346  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
347  rejectionReason |= kTimeCut;
348  return kFALSE;
349  }
350 
351  if (fExoticCut && clus->GetIsExotic()) {
352  rejectionReason |= kExoticCut;
353  return kFALSE;
354  }
355 
356  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
357  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
358  rejectionReason |= kEnergyCut;
359  return kFALSE;
360  }
361  }
362 
363  return kTRUE;
364 }
365 
371 {
372  UInt_t rejectionReason = 0;
373  Int_t nClus = 0;
374  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
375  AliVCluster *clust = this->GetCluster(iclust);
376  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
377  }
378  return nClus;
379 }
380 
387 {
388  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
389  return fUserDefEnergyCut[t];
390  }
391  else {
392  return fMinE;
393  }
394 }
395 
402 {
403  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
404  fUserDefEnergyCut[t] = cut;
405  }
406  else {
407  fMinE = cut;
408  }
409 }
410 
419 void AliClusterContainer::SetArray(const AliVEvent * event)
420 {
421  AliEmcalContainer::SetArray(event);
422 
423  // Register TClonesArray in index map
425 }
426 
433  return AliClusterIterableContainer(this, false);
434 }
435 
442  return AliClusterIterableContainer(this, true);
443 }
444 
451  return AliClusterIterableMomentumContainer(this, false);
452 }
453 
460  return AliClusterIterableMomentumContainer(this, true);
461 }
462 
463 const char* AliClusterContainer::GetTitle() const
464 {
465  static TString clusterString;
466 
468 
469  if (Ecut == 0) {
470  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
471  }
472  else if (Ecut < 1.0) {
473  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
474  }
475  else {
476  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
477  }
478 
479  return clusterString.Data();
480 }
481 
482 
483 /******************************************
484  * Unit tests *
485  ******************************************/
486 
487 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
488  std::vector<AliVCluster *> reference, variation;
489  AliVCluster *test = NULL;
490  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
491  test = cont->GetCluster(iclust);
492  if(!iteratorType){
493  UInt_t rejectionReason = 0;
494  if(!cont->AcceptCluster(test, rejectionReason)) continue;
495  }
496  reference.push_back(test);
497  }
498 
499  if(!iteratorType){
500  // test accept iterator
501  for(auto cluster : cont->accepted()){
502  variation.push_back(cluster);
503  }
504  } else {
505  // test all iterator
506  for(auto cluster : cont->all()){
507  variation.push_back(cluster);
508  }
509  }
510 
511  if(verbose){
512  // Printing cluster addresses:
513  if(reference.size() < 30){
514  std::cout << "Clusters in reference container: " << std::endl;
515  std::cout << "===========================================" << std::endl;
516  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
517  std::cout << "Address: " << *refit << std::endl;
518  }
519  }
520  if(variation.size() < 30){
521  std::cout << "Clusters in test container: " << std::endl;
522  std::cout << "===========================================" << std::endl;
523  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
524  std::cout << "Address: " << *varit << std::endl;
525  }
526  }
527  }
528 
529  int testresult = 0;
530  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
531  bool failure = false;
532  // first test: cleck if all clusters are found by the iterator
533  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
534  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
535  if(verbose)
536  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
537  failure = true;
538  break;
539  }
540  }
541  if(!failure){
542  // second test: check if there are no extra clusters found
543  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
544  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
545  if(verbose)
546  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
547  failure = true;
548  break;
549  }
550  }
551  if(failure) testresult = 2;
552  } else {
553  testresult = 1;
554  }
555 
556  if(verbose){
557  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
558  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
559  std::cout << "Test result: " << testresult << std::endl;
560  std::cout << "-----------------------------------------------------------------------" << std::endl;
561  }
562  return testresult;
563 }
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
void SetArray(const AliVEvent *event)
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()