AliPhysics  3337bb0 (3337bb0)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliParticleContainer.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 
23 #include "AliTLorentzVector.h"
24 #include "AliParticleContainer.h"
25 
27 ClassImp(AliParticleContainer);
29 
30 // Properly instantiate the object
32 
37  AliEmcalContainer(),
38  fMinDistanceTPCSectorEdge(-1),
39  fChargeCut(kNoChargeCut),
40  fGeneratorIndex(-1)
41 {
42  fBaseClassName = "AliVParticle";
43  SetClassName("AliVParticle");
44 }
45 
51  AliEmcalContainer(name),
52  fMinDistanceTPCSectorEdge(-1),
53  fChargeCut(kNoChargeCut),
54  fGeneratorIndex(-1)
55 {
56  fBaseClassName = "AliVParticle";
57  SetClassName("AliVParticle");
58 }
59 
66 AliVParticle* AliParticleContainer::GetLeadingParticle(const char* opt)
67 {
68  TString option(opt);
69  option.ToLower();
70 
71  Int_t tempID = fCurrentID;
72  ResetCurrentID();
73 
74  AliVParticle *partMax = GetNextAcceptParticle();
75  AliVParticle *part = 0;
76 
77  if (option.Contains("p")) {
78  while ((part = GetNextAcceptParticle())) {
79  if (part->P() > partMax->P()) partMax = part;
80  }
81  }
82  else {
83  while ((part = GetNextAcceptParticle())) {
84  if (part->Pt() > partMax->Pt()) partMax = part;
85  }
86  }
87 
88  fCurrentID = tempID;
89 
90  return partMax;
91 }
92 
99 {
100  //
101  if (i == -1) i = fCurrentID;
102  if (i < 0 || i >= this->fClArray->GetEntriesFast()) return 0;
103  AliVParticle *vp = static_cast<AliVParticle*>(fClArray->At(i));
104  return vp;
105 }
106 
114 {
115  UInt_t rejectionReason = 0;
116  if (i == -1) i = fCurrentID;
117  if (AcceptParticle(i, rejectionReason)) {
118  return GetParticle(i);
119  }
120  else {
121  AliDebug(2,"Particle not accepted.");
122  return 0;
123  }
124 }
125 
133 {
134  const Int_t n = GetNEntries();
135  AliVParticle *p = 0;
136  do {
137  fCurrentID++;
138  if (fCurrentID >= n) break;
139  p = GetAcceptParticle(fCurrentID);
140  } while (!p);
141 
142  return p;
143 }
144 
152 {
153  const Int_t n = GetNEntries();
154  AliVParticle *p = 0;
155  do {
156  fCurrentID++;
157  if (fCurrentID >= n) break;
158  p = GetParticle(fCurrentID);
159  } while (!p);
160 
161  return p;
162 }
163 
173 Bool_t AliParticleContainer::GetMomentumFromParticle(TLorentzVector &mom, const AliVParticle* part, Double_t mass) const
174 {
175  if (part && part->Eta() < 1e6 && part->Eta() > -1e6) { // protection against FPE in sinh(eta)
176  if (mass < 0) mass = part->M();
177  mom.SetPtEtaPhiM(part->Pt(), part->Eta(), part->Phi(), mass);
178  return kTRUE;
179  }
180  else {
181  mom.SetPtEtaPhiM(0, 0, 0, 0);
182  return kFALSE;
183  }
184 }
185 
193 Bool_t AliParticleContainer::GetMomentumFromParticle(TLorentzVector &mom, const AliVParticle* part) const
194 {
195  return GetMomentumFromParticle(mom,part,fMassHypothesis);
196 }
197 
207 Bool_t AliParticleContainer::GetMomentum(TLorentzVector &mom, Int_t i) const
208 {
209  if (i == -1) i = fCurrentID;
210  AliVParticle *vp = GetParticle(i);
211  return GetMomentumFromParticle(mom, vp);
212 }
213 
224 {
225  AliVParticle *vp = GetNextParticle();
226  return GetMomentumFromParticle(mom, vp);
227 }
228 
240 {
241  if (i == -1) i = fCurrentID;
242  AliVParticle *vp = GetAcceptParticle(i);
243  return GetMomentumFromParticle(mom, vp);
244 }
245 
256 {
257  AliVParticle *vp = GetNextAcceptParticle();
258  return GetMomentumFromParticle(mom, vp);
259 }
260 
271 Bool_t AliParticleContainer::AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
272 {
273  Bool_t r = ApplyParticleCuts(vp, rejectionReason);
274  if (!r) return kFALSE;
275 
276  AliTLorentzVector mom;
277 
278  Int_t id = fClArray->IndexOf(vp);
279  if (id >= 0) {
280  GetMomentum(mom, id);
281  }
282  else {
283  GetMomentumFromParticle(mom, vp);
284  }
285 
286  return ApplyKinematicCuts(mom, rejectionReason);
287 }
288 
300 {
301  Bool_t r = ApplyParticleCuts(GetParticle(i), rejectionReason);
302  if (!r) return kFALSE;
303 
304  AliTLorentzVector mom;
305  GetMomentum(mom, i);
306 
307  return ApplyKinematicCuts(mom, rejectionReason);
308 }
309 
322 Bool_t AliParticleContainer::ApplyParticleCuts(const AliVParticle* vp, UInt_t &rejectionReason) const
323 {
324  if (!vp) {
325  rejectionReason |= kNullObject;
326  return kFALSE;
327  }
328 
329  if (vp->TestBits(fBitMap) != (Int_t)fBitMap) {
330  rejectionReason |= kBitMapCut;
331  return kFALSE;
332  }
333 
334  if (fMinMCLabel >= 0 && TMath::Abs(vp->GetLabel()) < fMinMCLabel) {
335  rejectionReason |= kMCLabelCut;
336  return kFALSE;
337  }
338 
339  if (fMaxMCLabel >= 0 && TMath::Abs(vp->GetLabel()) > fMaxMCLabel) {
340  rejectionReason |= kMCLabelCut;
341  return kFALSE;
342  }
343 
344  switch (fChargeCut) {
345  case kCharged:
346  if (vp->Charge() == 0) {
347  rejectionReason |= kChargeCut;
348  return kFALSE;
349  }
350  break;
351 
352  case kNeutral:
353  if (vp->Charge() != 0) {
354  rejectionReason |= kChargeCut;
355  return kFALSE;
356  }
357  break;
358 
359  case kPositiveCharge:
360  if (vp->Charge() <= 0) {
361  rejectionReason |= kChargeCut;
362  return kFALSE;
363  }
364  break;
365 
366  case kNegativeCharge:
367  if (vp->Charge() >= 0) {
368  rejectionReason |= kChargeCut;
369  return kFALSE;
370  }
371  break;
372 
373  default:
374  break;
375  }
376 
377  if (fGeneratorIndex >= 0 && fGeneratorIndex != vp->GetGeneratorIndex()) {
378  rejectionReason |= kMCGeneratorCut;
379  return kFALSE;
380  }
381 
382  return kTRUE;
383 }
384 
397 {
399  const Double_t pi = TMath::Pi();
400  const Double_t kSector = pi/9;
401  Double_t phiDist = TMath::Abs(mom.Phi() - TMath::FloorNint(mom.Phi()/kSector)*kSector);
402  if(phiDist<fMinDistanceTPCSectorEdge) {
403  rejectionReason |= kMinDistanceTPCSectorEdgeCut;
404  return kFALSE;
405  }
406  }
407 
408  return AliEmcalContainer::ApplyKinematicCuts(mom, rejectionReason);
409 }
410 
418 {
419  Int_t nPart = 0;
420  for(int ipart = 0; ipart < this->GetNParticles(); ipart++){
421  UInt_t rejectionReason = 0;
422  if(this->AcceptParticle(ipart, rejectionReason)) nPart++;
423  }
424  return nPart;
425 }
426 
433 {
434  static TString trackString;
435 
436  if (GetMinPt() == 0) {
437  trackString = TString::Format("%s_pT0000", GetArrayName().Data());
438  }
439  else if (GetMinPt() < 1.0) {
440  trackString = TString::Format("%s_pT0%3.0f", GetArrayName().Data(), GetMinPt()*1000.0);
441  }
442  else {
443  trackString = TString::Format("%s_pT%4.0f", GetArrayName().Data(), GetMinPt()*1000.0);
444  }
445 
446  return trackString.Data();
447 }
448 
457 void AliParticleContainer::SetArray(const AliVEvent * event)
458 {
459  AliEmcalContainer::SetArray(event);
460 
461  // Register TClonesArray in index map
463 }
464 
471  return AliParticleIterableContainer(this, false);
472 }
473 
480  return AliParticleIterableContainer(this, true);
481 }
482 
489  return AliParticleIterableMomentumContainer(this, false);
490 }
491 
498  return AliParticleIterableMomentumContainer(this, true);
499 }
500 
501 /******************************************
502  * Unit tests *
503  ******************************************/
504 
505 int TestParticleContainerIterator(const AliParticleContainer *const cont, int iteratorType, bool verbose){
506  std::vector<AliVParticle *> reference, variation;
507  AliVParticle *test = NULL;
508  for(int iclust = 0; iclust < cont->GetNParticles(); iclust++){
509  test = cont->GetParticle(iclust);
510  if(!iteratorType){
511  UInt_t rejectionReason = 0;
512  if(!cont->AcceptParticle(test, rejectionReason)) continue;
513  }
514  reference.push_back(test);
515  }
516 
517  if(!iteratorType){
518  // test accept iterator
519  for(auto part : cont->accepted()){
520  variation.push_back(part);
521  }
522  } else {
523  // test all iterator
524  for(auto part : cont->all()){
525  variation.push_back(part);
526  }
527  }
528 
529  if(verbose) {
530  // Printing cluster addresses:
531  if(reference.size() < 30){
532  std::cout << "Paritcles in reference container: " << std::endl;
533  std::cout << "===========================================" << std::endl;
534  for(std::vector<AliVParticle *>::iterator refit = reference.begin(); refit != reference.end(); ++refit){
535  std::cout << "Address: " << *refit << std::endl;
536  }
537  }
538  if(variation.size() < 30){
539  std::cout << "Paritcles in test container: " << std::endl;
540  std::cout << "===========================================" << std::endl;
541  for(std::vector<AliVParticle *>::iterator varit = variation.begin(); varit != variation.end(); ++varit){
542  std::cout << "Address: " << *varit << std::endl;
543  }
544  }
545  }
546 
547  int testresult = 0;
548  // compare distributions: all particles in one vector needs to be found in the other vector and vice versa
549  bool failure = false;
550  // first test: cleck if all particles are found by the iterator
551  for(std::vector<AliVParticle *>::iterator clit = reference.begin(); clit != reference.end(); ++clit){
552  if(std::find(variation.begin(), variation.end(), *clit) == variation.end()) {
553  if(verbose)
554  std::cout << "Could not find particle with address " << *clit << " in test container" << std::endl;
555  failure = true;
556  break;
557  }
558  }
559  if(!failure){
560  // second test: check if there are no extra particles found
561  for(std::vector<AliVParticle *>::iterator clit = variation.begin(); clit != variation.end(); ++clit){
562  if(std::find(reference.begin(), reference.end(), *clit) == reference.end()) {
563  if(verbose)
564  std::cout << "Could not find particle with address " << *clit << " in reference container" << std::endl;
565  failure = true;
566  break;
567  }
568  }
569  if(failure) testresult = 2;
570  } else testresult = 1;
571  if(verbose) {
572  std::cout << "Unit test particle container, iterator type " << iteratorType << std::endl;
573  std::cout << "Number of expected particles: " << reference.size() << ", number of found particles: " << variation.size() << std::endl;
574  std::cout << "Test result: " << testresult << std::endl;
575  std::cout << "-----------------------------------------------------------------------" << std::endl;
576  }
577  return testresult;
578 }
const Double_t kSector
virtual AliVParticle * GetNextAcceptParticle()
double Double_t
Definition: External.C:58
void SetArray(const AliVEvent *event)
virtual Bool_t GetMomentum(TLorentzVector &mom, Int_t i) const
const char * GetTitle() const
const AliParticleIterableContainer all() const
Double_t mass
Int_t GetNParticles() const
Declaration of class AliTLorentzVector.
virtual Bool_t ApplyParticleCuts(const AliVParticle *vp, UInt_t &rejectionReason) const
EChargeCut_t fChargeCut
select particles according to their charge
virtual Bool_t GetNextMomentum(TLorentzVector &mom)
Container for particles within the EMCAL framework.
virtual Bool_t GetNextAcceptMomentum(TLorentzVector &mom)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
virtual AliVParticle * GetLeadingParticle(const char *opt="")
virtual AliVParticle * GetParticle(Int_t i=-1) const
EMCALIterableContainer::AliEmcalIterableContainerT< AliVParticle, EMCALIterableContainer::operator_star_object< AliVParticle > > AliParticleIterableContainer
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
virtual Bool_t AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
Double_t fMinDistanceTPCSectorEdge
require minimum distance to edge of TPC sector edge
static AliEmcalContainerIndexMap< TClonesArray, AliVParticle > fgEmcalContainerIndexMap
! Mapping from containers to indices
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)
Short_t fGeneratorIndex
select MC particles with generator index (default = -1 = switch off selection)
virtual Bool_t GetMomentumFromParticle(TLorentzVector &mom, const AliVParticle *part, Double_t mass) const
const AliParticleIterableMomentumContainer accepted_momentum() const
virtual Bool_t ApplyKinematicCuts(const AliTLorentzVector &mom, UInt_t &rejectionReason) const
const AliParticleIterableContainer accepted() const
void test(int runnumber=195345)
int TestParticleContainerIterator(const AliParticleContainer *const cont, int iteratorType, bool verbose)
EMCALIterableContainer::AliEmcalIterableContainerT< AliVParticle, EMCALIterableContainer::operator_star_pair< AliVParticle > > AliParticleIterableMomentumContainer
const Double_t pi
bool Bool_t
Definition: External.C:53
Int_t GetNAcceptedParticles() const
virtual Bool_t GetAcceptMomentum(TLorentzVector &mom, Int_t i) const
const AliParticleIterableMomentumContainer all_momentum() const
virtual AliVParticle * GetNextParticle()