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