JETAN/AliJetFinder.cxx

   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
  16 //---------------------------------------------------------------------
  17 // Jet finder base class
  18 // manages the search for jets
  19 // Author: jgcn@mda.cinvestav.mx
  20 //---------------------------------------------------------------------
  21
  22 #include <Riostream.h>
  23 #include <TFile.h>
  24 #include "AliGenPythiaEventHeader.h"
  25 #include "AliJetFinder.h"
  26 #include "AliJet.h"
  27 #include "AliJetReader.h"
  28 #include "AliJetReaderHeader.h"
  29 #include "AliJetControlPlots.h"
  30 #include "AliLeading.h"
  31 #include "AliHeader.h"
  32
  33 ClassImp(AliJetFinder)
  34
  35
  36 AliJetFinder::AliJetFinder():
  37   fPlotMode(kFALSE),
  38   fJets(0),
  39   fGenJets(0),
  40   fLeading(0),
  41   fReader(0x0),
  42   fPlots(0x0),
  43   fOut(0x0)
  44
  45 {
  46   // Constructor
  47   //fOut     = 0x0;
  48   fJets    = new AliJet();
  49   fGenJets = new AliJet();
  50   fLeading = new AliLeading();
  51   //fReader  = 0x0;
  52   //fPlots   = 0x0;
  53   //SetPlotMode(kFALSE);
  54 }
  55
  56 ////////////////////////////////////////////////////////////////////////
  57
  58 AliJetFinder::~AliJetFinder()
  59 {
  60   // destructor
  61   // here reset and delete jets
  62   fJets->ClearJets();
  63   delete fJets;
  64   fGenJets->ClearJets();
  65   delete fGenJets;
  66   // close file
  67   if (fOut) {
  68     fOut->Close();
  69     fOut->Delete();
  70   }
  71   delete fOut;
  72   // reset and delete control plots
  73   if (fPlots) delete fPlots;
  74 }
  75
  76 ////////////////////////////////////////////////////////////////////////
  77
  78 void AliJetFinder::SetOutputFile(const char *name)
  79 {
  80   // opens output file
  81   fOut = new TFile(name,"recreate");
  82 }
  83
  84 ////////////////////////////////////////////////////////////////////////
  85
  86 void AliJetFinder::PrintJets()
  87 {
  88   // Print jet information
  89   cout << " Jets found with jet algorithm:" << endl;
  90   fJets->PrintJets();
  91   cout << " Jets found by pythia:" << endl;
  92   fGenJets->PrintJets();
  93 }
  94
  95 ////////////////////////////////////////////////////////////////////////
  96
  97 void AliJetFinder::SetPlotMode(Bool_t b)
  98 {
  99   // Sets the plotting mode
 100   fPlotMode=b;
 101   if (b && !fPlots) fPlots = new AliJetControlPlots();
 102 }
 103
 104 ////////////////////////////////////////////////////////////////////////
 105
 106 void AliJetFinder::WriteJetsToFile(Int_t i)
 107 {
 108   // Writes the jets to file
 109   fOut->cd();
 110   char hname[30];
 111   sprintf(hname,"TreeJ%d",i);
 112   TTree* jetT = new TTree(hname,"AliJet");
 113   jetT->Branch("FoundJet",&fJets,1000);
 114   jetT->Branch("GenJet",&fGenJets,1000);
 115   jetT->Branch("LeadingPart",&fLeading,1000);
 116   jetT->Fill();
 117   jetT->Write(hname);
 118   delete jetT;
 119 }
 120
 121 ////////////////////////////////////////////////////////////////////////
 122
 123 void AliJetFinder::WriteRHeaderToFile()
 124 {
 125   // write reader header
 126   fOut->cd();
 127   AliJetReaderHeader *rh = fReader->GetReaderHeader();
 128   rh->Write();
 129 }
 130
 131 ////////////////////////////////////////////////////////////////////////
 132
 133 void AliJetFinder::GetGenJets()
 134 {
 135   // Get the generated jet information from mc header
 136   AliHeader* alih = fReader->GetAliHeader();
 137   if (alih == 0) return;
 138   AliGenEventHeader * genh = alih->GenEventHeader();
 139   if (genh == 0) return;
 140   Int_t nj =((AliGenPythiaEventHeader*)genh)->NTriggerJets();
 141   Int_t* m = new Int_t[nj];
 142   Int_t* k = new Int_t[nj];
 143   for (Int_t i=0; i< nj; i++) {
 144     Float_t p[4];
 145     ((AliGenPythiaEventHeader*)genh)->TriggerJet(i,p);
 146     fGenJets->AddJet(p[0],p[1],p[2],p[3]);
 147     m[i]=1;
 148     k[i]=i;
 149   }
 150   fGenJets->SetNinput(nj);
 151   fGenJets->SetMultiplicities(m);
 152   fGenJets->SetInJet(k);
 153 }
 154
 155 ////////////////////////////////////////////////////////////////////////
 156
 157 void AliJetFinder::Run()
 158 {
 159   // do some initialization
 160   Init();
 161
 162   // connect files
 163   fReader->OpenInputFiles();
 164
 165   // write headers
 166   if (fOut) {
 167       fOut->cd();
 168       WriteRHeaderToFile();
 169       WriteJHeaderToFile();
 170   }
 171   // loop over events
 172   Int_t nFirst,nLast;
 173   nFirst = fReader->GetReaderHeader()->GetFirstEvent();
 174   nLast = fReader->GetReaderHeader()->GetLastEvent();
 175   // loop over events
 176   for (Int_t i=nFirst;i<nLast;i++) {
 177       fReader->FillMomentumArray(i);
 178       fLeading->FindLeading(fReader);
 179       GetGenJets();
 180       FindJets();
 181       if (fOut) WriteJetsToFile(i);
 182       if (fPlots) fPlots->FillHistos(fJets);
 183       fLeading->Reset();
 184       fGenJets->ClearJets();
 185       Reset();
 186   }
 187   // write out
 188   if (fPlots) {
 189       fPlots->Normalize();
 190       fPlots->PlotHistos();
 191   }
 192   if (fOut) {
 193       fOut->cd();
 194       fPlots->Write();
 195       fOut->Close();
 196   }
 197 }