JETAN/AliTkJetTrigger.cxx

   1 // $Id$
   2
   3 /****************************************************************************
   4  * AliTkJetTrigger.cxx                                                         *
   5  * Test version of jet-trigger, includes many analysis functions to tune    *
   6  * paramters                                                                *
   7  * Thorsten Kollegger <kollegge@ikf.physik.uni-frankfurt.de>                *
   8  ***************************************************************************/
   9
  10 #include <Riostream.h>
  11 #include <TFile.h>
  12 #include <TROOT.h>
  13 #include <TParticle.h>
  14 #include <TClonesArray.h>
  15 #include <TTree.h>
  16 #include "AliTkChargedJetFinder.h"
  17 #include "AliTkJetTriggerDecision.h"
  18 #include "AliTkJetTriggerEvent.h"
  19 #include "AliTkJetTrigger.h"
  20
  21 void AliTkJetTrigger::init() {
  22   // open file and create tree...
  23   if (getEvOutFilename() != NULL) {
  24     cout << "Writing events to " << getEvOutFilename() << endl;
  25     mEvOutFile = new TFile(getEvOutFilename(),"RECREATE");
  26     mEvOutEvent = new AliTkJetTriggerEvent();
  27     mEvOutTree = new TTree("jetTrigger","");
  28     mEvOutTree->Branch("event","AliTkJetTriggerEvent",&mEvOutEvent,32000,99);
  29   }
  30 }
  31
  32 void AliTkJetTrigger::initEvent() {
  33   // can not be implemented, we need the particles...
  34 }
  35
  36 void AliTkJetTrigger::initEvent(TClonesArray *particles, Int_t type) {
  37   setParticleType(type);
  38   initEvent(particles);
  39 }
  40
  41 void AliTkJetTrigger::initEvent(TClonesArray *particles) {
  42   const UInt_t mParticlesDefaultSize = getExpectedParticleNumber();
  43   // make a copy of all accepted particles
  44   // just to be sure that they are not modified...
  45
  46   clearParticles();
  47   if (!particles) {
  48     return;
  49   }
  50   mParticles = new TClonesArray("TParticle",mParticlesDefaultSize);
  51   TIterator *iter = particles->MakeIterator();
  52   TParticle *particle;
  53   UInt_t pos = 0;
  54   while ((particle = (TParticle *) iter->Next()) != NULL) {
  55     if (isParticleAccepted(particle)) {
  56       if(particle->Pt()>fParticleMinPt)
  57         new ((*mParticles)[pos++]) TParticle(*particle);
  58     }
  59   }
  60   delete iter;
  61   createSeedPoints();
  62 }
  63
  64 void AliTkJetTrigger::make() {
  65   // fill the matrix containing trigger decision for all different cone radia,
  66   // nParticles and PtThresholds...
  67   // loop over all seedpoints...
  68   cout << "generated " << seedPoints.size() << " seedpoints" << endl;
  69   for (list<AliTkEtaPhiVector>::const_iterator seedIter = seedPoints.begin();
  70        seedIter != seedPoints.end(); ++seedIter) {
  71     AliTkEtaPhiVector coneCenter(*seedIter);
  72     for (list<Float_t>::const_iterator radiusIter = mConeRadia.begin();
  73          radiusIter != mConeRadia.end(); ++radiusIter) {
  74       // get particles in cone
  75       list<TParticle*> *myParticles = getParticlesInCone(coneCenter,
  76                                                          *radiusIter);
  77       // loop over all PtThresholds...
  78       for (list<Float_t>::const_iterator PtThreshold =mConePtThreshold.begin();
  79            PtThreshold != mConePtThreshold.end(); ++PtThreshold) {
  80         // loop over all nParticles thresholds...
  81         for (list<UInt_t>::const_iterator nParticles = mConeNParticles.begin();
  82              nParticles != mConeNParticles.end(); ++nParticles) {
  83           // what's the decision???
  84           if (decide(myParticles,*PtThreshold,*nParticles)) {
  85             AliTkJetTriggerDecision *d = new AliTkJetTriggerDecision();
  86             d->setConeRadius(*radiusIter);
  87             d->setPtThreshold(*PtThreshold);
  88             d->setNParticles(*nParticles);
  89             mEvOutEvent->addDecision(d);
  90             // new change, requires testing
  91             delete d;
  92           }
  93         } // end of nParticle threshold loop
  94       } // end of PtThreshold loop
  95       myParticles->erase(myParticles->begin(),myParticles->end());
  96       delete myParticles;
  97     } // end of cone radia loop
  98   } // end of seed point loop
  99 }
 100
 101 void AliTkJetTrigger::finishEvent() {
 102   // let's write out all triggered settings...
 103   // not finish yet...
 104   clearParticles();
 105   if (mEvOutFile) {
 106     mEvOutTree->Fill();
 107   }
 108   mEvOutEvent->clear();
 109   seedPoints.erase(seedPoints.begin(),seedPoints.end());
 110 }
 111
 112 void AliTkJetTrigger::finish() {
 113   // write object containing settings
 114   // close file
 115   if (mEvOutFile) {
 116     if (!mEvOutFile->IsOpen()) { cerr << "WHAT?" << endl; }
 117     mEvOutFile->cd();
 118     mEvOutTree->Write();
 119     mEvOutFile->Close();
 120   }
 121 }
 122
 123 void AliTkJetTrigger::setDefaults() {
 124   fParticleMinPt=0;
 125   fParticleMinPtSeed=1;
 126   setParticleType(3);
 127   addConeRadius(0.7);
 128   addConeRadius(0.5);
 129   addConeRadius(0.3);
 130   addConeRadius(0.2);
 131   addConeRadius(0.1);
 132   addConeNParticles(1);
 133   addConeNParticles(2);
 134   addConeNParticles(3);
 135   addConeNParticles(4);
 136   addConeNParticles(5);
 137   addConePtThreshold(1.);
 138   addConePtThreshold(2.);
 139   addConePtThreshold(3.);
 140   addConePtThreshold(4.);
 141   addConePtThreshold(5.);
 142
 143   setEvOutFilename("$JF_DATADIR/trigger.evout.root");
 144   setExpectedParticleNumber(30000);
 145 }
 146
 147 void AliTkJetTrigger::setParticleType(Int_t type) {
 148   mParticleType = type;
 149   // current used particles types:
 150   // 1 == TParticle from TPythia6
 151   // 2 == TParticle from THijing
 152   // 3 == TParticle from TPythia6 + ALICE acceptance approximation
 153   // 4 == TParticle from THijing  + ALICE acceptance approximation
 154 }
 155
 156 void AliTkJetTrigger::addConeRadius(Float_t radius) {
 157   mConeRadia.push_back(radius);
 158 }
 159
 160 void AliTkJetTrigger::addConeNParticles(UInt_t nParticles) {
 161   mConeNParticles.push_back(nParticles);
 162 }
 163
 164 void AliTkJetTrigger::addConePtThreshold(Float_t ptThreshold) {
 165   mConePtThreshold.push_back(ptThreshold);
 166 }
 167
 168 void AliTkJetTrigger::mixParticles(TClonesArray *particles, Int_t type) {
 169   // mixes particles to the existing event
 170   if (!mParticles) {
 171     return;
 172   }
 173   if (!particles) {
 174     return;
 175   }
 176   TParticle *particle = NULL;
 177   UInt_t pos = 0;
 178   pos = mParticles->GetLast()+1;
 179   TIterator *iter = particles->MakeIterator();
 180   while ((particle = (TParticle *) iter->Next()) != NULL) {
 181     if (isParticleAccepted(particle,type)) {
 182       new ((*mParticles)[pos++]) TParticle(*particle);
 183     }
 184   }
 185   delete iter;
 186   // we need to generate new seedPoints...
 187   seedPoints.erase(seedPoints.begin(),seedPoints.end());
 188   createSeedPoints();
 189 }
 190
 191 void AliTkJetTrigger::createSeedPoints() {
 192   TParticle *particle;
 193   TIterator *iter = mParticles->MakeIterator();
 194   while ((particle = (TParticle *)iter->Next()) != NULL) {
 195     // particles are checked for acceptance during initEvent, mixParticles
 196     //    if (isParticleAccepted(particle)) {
 197       if (particle->Pt() > fParticleMinPtSeed) {
 198         AliTkEtaPhiVector v(particle->Eta(),particle->Phi());
 199         seedPoints.push_back(v);
 200       }
 201       //    }
 202   }
 203   delete iter;
 204 }
 205
 206 Bool_t AliTkJetTrigger::isParticleAccepted(TParticle *particle) {
 207   return isParticleAccepted(particle,getParticleType());
 208 }
 209
 210 Bool_t AliTkJetTrigger::isParticleAccepted(TParticle *particle, Int_t type) {
 211   switch (type) {
 212   case 1:
 213     return isParticleAcceptedPythia(particle);
 214     break;
 215   case 2:
 216     return isParticleAcceptedHijing(particle);
 217     break;
 218   case 3:
 219     return isParticleAcceptedPythiaAliceGeo(particle);
 220     break;
 221   case 4:
 222     return isParticleAcceptedHijingAliceGeo(particle);
 223     break;
 224   };
 225   return kFALSE;
 226 }
 227
 228 Bool_t AliTkJetTrigger::isParticleAcceptedPythia(TParticle *particle) {
 229   // check for stable particle
 230   if (particle->GetStatusCode() != 1) {
 231     return kFALSE;
 232   }
 233   return kTRUE;
 234 }
 235
 236 Bool_t AliTkJetTrigger::isParticleAcceptedPythiaAliceGeo(TParticle *particle) {
 237   // check if valid PYTHIA particle
 238   if (!isParticleAcceptedPythia(particle)) {
 239     return kFALSE;
 240   }
 241   // check if particle is in |eta| < 1
 242   if (TMath::Abs(particle->Eta()) > 1) {
 243     return kFALSE;
 244   }
 245   // check if charged particle
 246   TParticlePDG *partPDG = particle->GetPDG();
 247   if (partPDG->Charge() == 0) {
 248     return kFALSE;
 249   }
 250   return kTRUE;
 251 }
 252
 253 Bool_t AliTkJetTrigger::isParticleAcceptedHijing(TParticle */*particle*/) {
 254   return kTRUE;
 255 }
 256
 257 Bool_t AliTkJetTrigger::isParticleAcceptedHijingAliceGeo(TParticle *particle) {
 258   if (!isParticleAcceptedHijing(particle)) {
 259     return kFALSE;
 260   }
 261   // check if particle is in |eta| < 1
 262   if (TMath::Abs(particle->Eta()) > 1) {
 263     return kFALSE;
 264   }
 265   // check if charged particle
 266   TParticlePDG *partPDG = particle->GetPDG();
 267   if (partPDG->Charge() == 0) {
 268     return kFALSE;
 269   }
 270   return kTRUE;
 271 }
 272
 273 void AliTkJetTrigger::clearParticles() {
 274   // delete my particle copy...
 275   if (mParticles) {
 276     mParticles->Delete();
 277     delete mParticles;
 278   }
 279   mParticles = NULL;
 280 }
 281
 282 list<TParticle*> *AliTkJetTrigger::getParticlesInCone(AliTkEtaPhiVector center,
 283                                                    Float_t radius) {
 284   return getParticlesInCone(center,radius,mParticles);
 285 }
 286
 287 list<TParticle*> *AliTkJetTrigger::getParticlesInCone(AliTkEtaPhiVector center,
 288                                                    Float_t radius,
 289                                                    TClonesArray *particles) {
 290   // returns a list with all particles which are in a cone with radius around
 291   // center
 292   // one can easily speed up this part by saving the last cone+parameters and
 293   // calculate the new one from this one if there is some overlap
 294   list<TParticle*> *myParticles = new list<TParticle*>;
 295   Float_t radiusSq = radius*radius;
 296   TIterator *iter = particles->MakeIterator();
 297   TParticle *particle;
 298   while ((particle = (TParticle *)iter->Next()) != NULL) {
 299     AliTkEtaPhiVector pos(particle->Eta(),particle->Phi());
 300     if (center.diffSq(pos) < radiusSq) {
 301       myParticles->push_back(particle);
 302     }
 303   }
 304   delete iter;
 305   return myParticles;
 306 }
 307
 308 Bool_t AliTkJetTrigger::decide(list<TParticle*> *particles, Float_t ptThreshold,
 309                             UInt_t nParticles) {
 310   UInt_t n = 0;
 311   for (list<TParticle*>::const_iterator iter = particles->begin();
 312        iter != particles->end(); ++iter) {
 313     if ((*iter)->Pt() > ptThreshold) {
 314       n++;
 315     }
 316   }
 317   if (n >= nParticles) {
 318     return kTRUE;
 319   }
 320   return kFALSE;
 321 }
 322
 323 void AliTkJetTrigger::setEvOutFilename(Char_t *filename) {
 324   if (!mEvOutName) {
 325     mEvOutName = new Char_t[4096];
 326   }
 327   strcpy(mEvOutName,filename);
 328 }
 329
 330
 331 ClassImp(AliTkJetTrigger)