JETAN/AliJFKtJetFinder.cxx

   1 // $Id$
   2
   3 /* Kt Jet finder class ala hep-ex/0005012
   4    for CDF/D0 Run II */
   5
   6 #include <Riostream.h>
   7 #include <vector>
   8 #include <set>
   9
  10 #include <TClonesArray.h>
  11 #include <TIterator.h>
  12 #include <TParticle.h>
  13 #include <TParticlePDG.h>
  14 #include <TObjArray.h>
  15
  16 #include "AliJFPreCluster.h"
  17 #include "AliJFCluster.h"
  18 #include "AliJFClusterDifference.h"
  19 #include "AliJFKtJet.h"
  20 #include "AliJFKtJetFinder.h"
  21
  22 ClassImp(AliJFKtJetFinder)
  23
  24 AliJFKtJetFinder::AliJFKtJetFinder(Int_t n) : AliJFJetFinder(n)
  25 {
  26 }
  27
  28 AliJFKtJetFinder::~AliJFKtJetFinder()
  29 {
  30   Clean();
  31 }
  32
  33 Bool_t AliJFKtJetFinder::IsAcceptedParticle(TParticle *p)
  34 {
  35   if(p->GetStatusCode()%100!=1) return kFALSE;
  36
  37   Int_t pcode=p->GetPdgCode();
  38
  39   if((!fEM) && ((pcode==11)||(pcode==-11)||(pcode==22))) return kFALSE;
  40
  41   TParticlePDG *pdg=p->GetPDG();
  42   Float_t ch=pdg->Charge();
  43   if((!fCharged)&&(ch)) return kFALSE;
  44   if((!fNeutral)&&(!ch)) return kFALSE;
  45
  46   Float_t eta=p->Eta();
  47   if((eta<fEtaMin)||(eta>fEtaMax)) return kFALSE;
  48
  49   Float_t phi=p->Phi();
  50   if((phi<fPhiMin)||(phi>fPhiMax)) return kFALSE;
  51
  52   Float_t pt=p->Pt();
  53   if((pt<fPtMin)||(pt>fPtMax)) return kFALSE;
  54
  55   return kTRUE;
  56 }
  57
  58 Int_t AliJFKtJetFinder::Init(TClonesArray *particles)
  59 { //create precluster according to particles
  60   //each precluster corresponds to a particle
  61   //in the accepted range
  62   if(particles==NULL) return -1;
  63
  64   TIterator *iter=particles->MakeIterator();
  65   TParticle *p;
  66   Int_t ret=0;
  67
  68   while((p=(TParticle*)iter->Next()) != NULL){
  69     if(IsAcceptedParticle(p)){
  70       ret++;
  71
  72       //fPreClusterList.push_back(new AliJFPreCluster(p));
  73       fPreClusterList.push_back(new AliJFPreCluster(p->Px(),p->Py(),p->Pz(),-1,p));
  74     }
  75   }
  76
  77 #if 0
  78   for(vector<AliJFPreCluster*>::iterator c=fPreClusterList.begin();c!=fPreClusterList.end();c++){
  79     cout << *(*c) << endl;
  80   }
  81   //exit(1);
  82 #endif
  83
  84   for(vector<AliJFPreCluster*>::iterator i=fPreClusterList.begin();i!=fPreClusterList.end();i++){
  85     fClusterList.push_back(new AliJFCluster(*(*i)));
  86   }
  87
  88 #if 0
  89   for(vector<AliJFCluster*>::iterator c=fClusterList.begin();c!=fClusterList.end();c++){
  90     cout << *(*c) << endl;
  91   }
  92   //exit(1);
  93 #endif
  94
  95   AliJFClusterDifference diff;
  96   for(vector<AliJFCluster*>::iterator i=fClusterList.begin();i!=fClusterList.end();i++){
  97     for(vector<AliJFCluster*>::iterator j=i;j!=fClusterList.end();j++){
  98       diff.SetValues(*i,*j);
  99       fClusterDiffSet.insert(diff);
 100     }
 101   }
 102
 103 #if 0
 104   for(multiset<AliJFClusterDifference>::iterator pos=fClusterDiffSet.begin();pos!=fClusterDiffSet.end();pos++){
 105     cout << *pos << endl;
 106   }
 107   //exit(1);
 108 #endif
 109
 110   return ret;
 111 }
 112
 113 Int_t AliJFKtJetFinder::Run()
 114 {
 115   //loop over stored diffence objects until set is empty
 116   while(!fClusterDiffSet.empty()){
 117     //get first element of sorted set
 118     multiset<AliJFClusterDifference>::iterator pos=fClusterDiffSet.begin();
 119
 120     //cout << pos->GetDij() << endl;
 121     if((!pos->IsValidPointer())||(!pos->IsValidEntry())){ //delete old value left in set
 122       fClusterDiffSet.erase(pos);
 123
 124       continue;
 125     }
 126
 127     //get information on stored difference
 128     AliJFCluster *i=pos->GetI();
 129     AliJFCluster *j=pos->GetJ();
 130     bool isjetfound=pos->IsDiagonal();
 131     fClusterDiffSet.erase(pos); //take it out
 132
 133     if(isjetfound){ //found jet
 134       i->MarkIsJet();
 135
 136       if(fNJets==fNJetsMax){
 137         fNJetsMax+=fNJets;
 138         fJets.Expand(fNJetsMax);
 139       }
 140
 141       fJets.AddAt(new AliJFKtJet(),fNJets);
 142       fJet=(AliJFKtJet*)fJets.At(fNJets);
 143       fJet->SetNJet(++fNJets);
 144       fJet->AddJet(i);
 145       fJet->Update();
 146
 147     } else { //combine cluster and make new difference objects
 148       i->CombineCluster(*j);
 149
 150       AliJFClusterDifference diff; //create and insert new difference objects
 151       for(vector<AliJFCluster*>::iterator c=fClusterList.begin();c!=fClusterList.end();c++){
 152         if((*c)->IsValid()) {
 153           diff.SetValues(i,*c);
 154           fClusterDiffSet.insert(diff);
 155         }
 156       }
 157     } //end if combine cluster
 158   } //end while
 159
 160   fJets.Sort();
 161   //fJets.Print();
 162
 163   return fNJets;
 164 }
 165
 166 void AliJFKtJetFinder::Debug()
 167 {
 168   for(vector<AliJFPreCluster*>::iterator c=fPreClusterList.begin();c!=fPreClusterList.end();c++){
 169     cout << *(*c) << endl;
 170   }
 171
 172   for(vector<AliJFCluster*>::iterator c=fClusterList.begin();c!=fClusterList.end();c++){
 173     //(*c)->Print();
 174   }
 175
 176 }
 177
 178 void AliJFKtJetFinder::Clean()
 179 {
 180   if(!fPreClusterList.empty()){
 181     for(vector<AliJFPreCluster*>::iterator i=fPreClusterList.begin();i!=fPreClusterList.end();i++){
 182       delete (*i);
 183     }
 184     fPreClusterList.erase(fPreClusterList.begin(),fPreClusterList.end());
 185   }
 186   if(!fClusterList.empty()){
 187     for(vector<AliJFCluster*>::iterator i=fClusterList.begin();i!=fClusterList.end();i++){
 188       delete (*i);
 189     }
 190     fClusterList.erase(fClusterList.begin(),fClusterList.end());
 191   }
 192
 193   if(!fClusterDiffSet.empty()) fClusterDiffSet.clear();
 194
 195   AliJFJetFinder::Clean();
 196 }