AliPhysics  6f1d526 (6f1d526)
 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 
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  fPhosMinNcells(0),
46  fPhosMinM02(0)
47 {
48  fBaseClassName = "AliVCluster";
49  SetClassName("AliVCluster");
50 
51  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
52  fUserDefEnergyCut[i] = 0.;
53  }
54 }
55 
61  AliEmcalContainer(name),
62  fClusTimeCutLow(-10),
63  fClusTimeCutUp(10),
64  fExoticCut(kTRUE),
65  fDefaultClusterEnergy(-1),
66  fIncludePHOS(kFALSE),
67  fPhosMinNcells(0),
68  fPhosMinM02(0)
69 {
70  fBaseClassName = "AliVCluster";
71  SetClassName("AliVCluster");
72 
73  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
74  fUserDefEnergyCut[i] = 0.;
75  }
76 }
77 
83 AliVCluster* AliClusterContainer::GetLeadingCluster(const char* opt)
84 {
85  TString option(opt);
86  option.ToLower();
87 
88  double (AliTLorentzVector::*momentum)() const = 0;
89 
90  if (option.Contains("e")) {
91  momentum = &AliTLorentzVector::E;
92  }
93  else {
94  momentum = &AliTLorentzVector::Et;
95  }
96 
98 
99  for (auto cluster : accepted_momentum()) {
100  if ((clusterMax.first.*momentum)() < (cluster.first.*momentum)()) {
101  clusterMax = cluster;
102  }
103  }
104 
105  return clusterMax.second;
106 }
107 
114 {
115  if(i<0 || i>fClArray->GetEntriesFast()) return 0;
116  AliVCluster *vp = static_cast<AliVCluster*>(fClArray->At(i));
117  return vp;
118 
119 }
120 
127 {
128  AliVCluster *vc = GetCluster(i);
129  if (!vc) return 0;
130 
131  UInt_t rejectionReason = 0;
132  if (AcceptCluster(vc, rejectionReason))
133  return vc;
134  else {
135  AliDebug(2,"Cluster not accepted.");
136  return 0;
137  }
138 }
139 
146 {
147  Int_t i = GetIndexFromLabel(lab);
148  return GetCluster(i);
149 }
150 
157 {
158  Int_t i = GetIndexFromLabel(lab);
159  return GetAcceptCluster(i);
160 }
161 
167 {
168  const Int_t n = GetNEntries();
169  AliVCluster *c = 0;
170  do {
171  fCurrentID++;
172  if (fCurrentID >= n) break;
173  c = GetAcceptCluster(fCurrentID);
174  } while (!c);
175 
176  return c;
177 }
178 
184 {
185  const Int_t n = GetNEntries();
186  AliVCluster *c = 0;
187  do {
188  fCurrentID++;
189  if (fCurrentID >= n) break;
190  c = GetCluster(fCurrentID);
191  } while (!c);
192 
193  return c;
194 }
195 
196 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc, Double_t mass) const
197 {
198  if (mass < 0) mass = 0;
199 
200  Double_t energy = 0;
201 
202  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
203  energy = vc->GetUserDefEnergy((AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
204  }
205  else {
206  energy = vc->E();
207  }
208 
209  Double_t p = TMath::Sqrt(energy*energy - mass*mass);
210 
211  Float_t pos[3];
212  vc->GetPosition(pos);
213 
214  pos[0]-=fVertex[0];
215  pos[1]-=fVertex[1];
216  pos[2]-=fVertex[2];
217 
218  Double_t r = TMath::Sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]) ;
219 
220  if (r > 1e-12) {
221  mom.SetPxPyPzE( p*pos[0]/r, p*pos[1]/r, p*pos[2]/r, energy) ;
222  }
223  else {
224  AliInfo("Null cluster radius, momentum calculation not possible");
225  return kFALSE;
226  }
227 
228  return kTRUE;
229 }
230 
231 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, const AliVCluster* vc) const
232 {
233  if (fMassHypothesis > 0) return GetMomentum(mom, vc, fMassHypothesis);
234 
235  if (vc) {
236  if (fDefaultClusterEnergy >= 0 && fDefaultClusterEnergy <= AliVCluster::kLastUserDefEnergy) {
237  vc->GetMomentum(mom, fVertex, (AliVCluster::VCluUserDefEnergy_t)fDefaultClusterEnergy);
238  }
239  else {
240  vc->GetMomentum(mom, fVertex);
241  }
242  }
243  else {
244  mom.SetPtEtaPhiM(0, 0, 0, 0.139);
245  return kFALSE;
246  }
247  return kTRUE;
248 }
249 
256 Bool_t AliClusterContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
257 {
258  AliVCluster *vc = GetCluster(i);
259  return GetMomentum(mom, vc);
260 }
261 
268 {
269  AliVCluster *vc = GetNextCluster();
270  return GetMomentum(mom, vc);
271 }
272 
280 {
281  AliVCluster *vc = GetAcceptCluster(i);
282  return GetMomentum(mom, vc);
283 }
284 
291 {
292  AliVCluster *vc = GetNextAcceptCluster();
293  return GetMomentum(mom, vc);
294 }
295 
297 {
298  Bool_t r = ApplyClusterCuts(GetCluster(i), rejectionReason);
299  if (!r) return kFALSE;
300 
301  AliTLorentzVector mom;
302  GetMomentum(mom, i);
303 
304  return ApplyKinematicCuts(mom, rejectionReason);
305 }
306 
307 Bool_t AliClusterContainer::AcceptCluster(const AliVCluster* clus, UInt_t &rejectionReason) const
308 {
309  Bool_t r = ApplyClusterCuts(clus, rejectionReason);
310  if (!r) return kFALSE;
311 
312  AliTLorentzVector mom;
313  GetMomentum(mom, clus);
314 
315  return ApplyKinematicCuts(mom, rejectionReason);
316 }
317 
323 Bool_t AliClusterContainer::ApplyClusterCuts(const AliVCluster* clus, UInt_t &rejectionReason) const
324 {
325  if (!clus) {
326  rejectionReason |= kNullObject;
327  return kFALSE;
328  }
329 
330  Bool_t bInAcceptance = clus->IsEMCAL();
331  if (fIncludePHOS) bInAcceptance = clus->IsEMCAL() || (clus->GetType() == AliVCluster::kPHOSNeutral);
332  if (!bInAcceptance) {
333  rejectionReason |= kIsEMCalCut;
334  return kFALSE;
335  }
336 
337  if (clus->TestBits(fBitMap) != (Int_t)fBitMap) {
338  rejectionReason |= kBitMapCut;
339  return kFALSE;
340  }
341 
342  if (fMinMCLabel >= 0 && TMath::Abs(clus->GetLabel()) < fMinMCLabel) {
343  rejectionReason |= kMCLabelCut;
344  return kFALSE;
345  }
346 
347  if (fMaxMCLabel >= 0 && TMath::Abs(clus->GetLabel()) > fMaxMCLabel) {
348  rejectionReason |= kMCLabelCut;
349  return kFALSE;
350  }
351 
352  if (clus->GetTOF() > fClusTimeCutUp || clus->GetTOF() < fClusTimeCutLow) {
353  rejectionReason |= kTimeCut;
354  return kFALSE;
355  }
356 
357  if (fExoticCut && clus->GetIsExotic()) {
358  rejectionReason |= kExoticCut;
359  return kFALSE;
360  }
361 
362  for (Int_t i = 0; i <= AliVCluster::kLastUserDefEnergy; i++) {
363  if (clus->GetUserDefEnergy((VCluUserDefEnergy_t)i) < fUserDefEnergyCut[i]) {
364  rejectionReason |= kEnergyCut;
365  return kFALSE;
366  }
367  }
368 
369  if (fIncludePHOS && clus->GetType() == AliVCluster::kPHOSNeutral) {
370  if (clus->GetNCells() < fPhosMinNcells) {
371  rejectionReason |= kExoticCut;
372  return kFALSE;
373  }
374 
375  if (clus->GetM02() < fPhosMinM02) {
376  rejectionReason |= kExoticCut;
377  return kFALSE;
378  }
379  }
380 
381  return kTRUE;
382 }
383 
389 {
390  UInt_t rejectionReason = 0;
391  Int_t nClus = 0;
392  for(int iclust = 0; iclust < this->fClArray->GetEntries(); ++iclust){
393  AliVCluster *clust = this->GetCluster(iclust);
394  if(this->AcceptCluster(clust, rejectionReason)) nClus++;
395  }
396  return nClus;
397 }
398 
405 {
406  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
407  return fUserDefEnergyCut[t];
408  }
409  else {
410  return fMinE;
411  }
412 }
413 
420 {
421  if (t >= 0 && t <= AliVCluster::kLastUserDefEnergy){
422  fUserDefEnergyCut[t] = cut;
423  }
424  else {
425  fMinE = cut;
426  }
427 }
428 
437 void AliClusterContainer::SetArray(const AliVEvent * event)
438 {
439  AliEmcalContainer::SetArray(event);
440 
441  // Register TClonesArray in index map
443 }
444 
451  return AliClusterIterableContainer(this, false);
452 }
453 
460  return AliClusterIterableContainer(this, true);
461 }
462 
469  return AliClusterIterableMomentumContainer(this, false);
470 }
471 
478  return AliClusterIterableMomentumContainer(this, true);
479 }
480 
481 const char* AliClusterContainer::GetTitle() const
482 {
483  static TString clusterString;
484 
486 
487  if (Ecut == 0) {
488  clusterString = TString::Format("%s_E0000", GetArrayName().Data());
489  }
490  else if (Ecut < 1.0) {
491  clusterString = TString::Format("%s_E0%3.0f", GetArrayName().Data(), Ecut*1000.0);
492  }
493  else {
494  clusterString = TString::Format("%s_E%4.0f", GetArrayName().Data(), Ecut*1000.0);
495  }
496 
497  return clusterString.Data();
498 }
499 
500 TString AliClusterContainer::GetDefaultArrayName(const AliVEvent * const ev) const {
501  if(ev->IsA() == AliAODEvent::Class()) return "caloClusters";
502  else if(ev->IsA() == AliESDEvent::Class()) return "CaloClusters";
503  else return "";
504 }
505 
506 
507 /******************************************
508  * Unit tests *
509  ******************************************/
510 
511 int TestClusterContainerIterator(const AliClusterContainer *const cont, int iteratorType, bool verbose){
512  std::vector<AliVCluster *> reference, variation;
513  AliVCluster *test = NULL;
514  for(int iclust = 0; iclust < cont->GetNClusters(); iclust++){
515  test = cont->GetCluster(iclust);
516  if(!iteratorType){
517  UInt_t rejectionReason = 0;
518  if(!cont->AcceptCluster(test, rejectionReason)) continue;
519  }
520  reference.push_back(test);
521  }
522 
523  if(!iteratorType){
524  // test accept iterator
525  for(auto cluster : cont->accepted()){
526  variation.push_back(cluster);
527  }
528  } else {
529  // test all iterator
530  for(auto cluster : cont->all()){
531  variation.push_back(cluster);
532  }
533  }
534 
535  if(verbose){
536  // Printing cluster addresses:
537  if(reference.size() < 30){
538  std::cout << "Clusters in reference container: " << std::endl;
539  std::cout << "===========================================" << std::endl;
540  for(std::vector<AliVCluster *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
541  std::cout << "Address: " << *refit << std::endl;
542  }
543  }
544  if(variation.size() < 30){
545  std::cout << "Clusters in test container: " << std::endl;
546  std::cout << "===========================================" << std::endl;
547  for(std::vector<AliVCluster *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
548  std::cout << "Address: " << *varit << std::endl;
549  }
550  }
551  }
552 
553  int testresult = 0;
554  // compare distributions: all clusters in one vector needs to be found in the other vector and vice versa
555  bool failure = false;
556  // first test: cleck if all clusters are found by the iterator
557  for(std::vector<AliVCluster *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
558  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
559  if(verbose)
560  std::cout << "Could not find cluster with address " << *clit << " in test container" << std::endl;
561  failure = true;
562  break;
563  }
564  }
565  if(!failure){
566  // second test: check if there are no extra clusters found
567  for(std::vector<AliVCluster *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
568  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
569  if(verbose)
570  std::cout << "Could not find cluster with address " << *clit << " in reference container" << std::endl;
571  failure = true;
572  break;
573  }
574  }
575  if(failure) testresult = 2;
576  } else {
577  testresult = 1;
578  }
579 
580  if(verbose){
581  std::cout << "Unit test cluster container, iterator type " << iteratorType << std::endl;
582  std::cout << "Number of expected clusters: " << reference.size() << ", number of found clusters: " << variation.size() << std::endl;
583  std::cout << "Test result: " << testresult << std::endl;
584  std::cout << "-----------------------------------------------------------------------" << std::endl;
585  }
586  return testresult;
587 }
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
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)
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)
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()
Int_t fPhosMinNcells
min number of phos cells per cluster