Updates needed for full jet reconstruction (charged + emcal) [Magali Estienne]
[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 /* $Id$ */
17
18 //---------------------------------------------------------------------
19 // Jet finder base class
20 // manages the search for jets 
21 // Authors: jgcn@mda.cinvestav.mx
22 //          andreas.morsch@cern.ch
23 //          magali.estienne@subatech.in2p3.fr
24 //---------------------------------------------------------------------
25
26 #include <Riostream.h>
27 #include <TFile.h>
28 #include <TClonesArray.h>
29 #include <TProcessID.h>
30
31 #include "AliJetFinder.h"
32 #include "AliJet.h"
33 #include "AliAODJet.h"
34 #include "AliJetReader.h"
35 #include "AliJetReaderHeader.h"
36 #include "AliJetControlPlots.h"
37 #include "AliLeading.h"
38 #include "AliAODEvent.h"
39 #include "AliJetUnitArray.h"
40
41 ClassImp(AliJetFinder)
42
43 AliJetFinder::AliJetFinder():
44     fTreeJ(0),
45     fPlotMode(kFALSE),
46     fJets(0),
47     fGenJets(0),
48     fLeading(0),
49     fReader(0x0),
50     fHeader(0x0),
51     fAODjets(0x0),
52     fNAODjets(0),
53     fPlots(0x0),
54     fOut(0x0)
55     
56 {
57   //
58   // Constructor
59   //
60
61   fJets    = new AliJet();
62   fGenJets = new AliJet();
63   fLeading = new AliLeading();
64   fAODjets = 0;
65 }
66
67 ////////////////////////////////////////////////////////////////////////
68 AliJetFinder::~AliJetFinder()
69 {
70   //
71   // Destructor
72   //
73
74   // Reset and delete jets
75   fJets->ClearJets();
76   delete fJets;
77   fGenJets->ClearJets();
78   delete fGenJets;
79   // close file
80   if (fOut) {
81     fOut->Close();
82     fOut->Delete();
83   }
84   delete fOut;
85   // reset and delete control plots
86   if (fPlots) delete fPlots;
87 }
88
89 ////////////////////////////////////////////////////////////////////////
90 void AliJetFinder::SetOutputFile(const char */*name*/)
91 {
92   //  opens output file 
93   //  fOut = new TFile(name,"recreate");
94 }
95
96 ////////////////////////////////////////////////////////////////////////
97 void AliJetFinder::PrintJets()
98 {
99   // Print jet information
100   cout << " Jets found with jet algorithm:" << endl;
101   fJets->PrintJets();
102   cout << " Jets found by pythia:" << endl;
103   fGenJets->PrintJets();
104 }
105
106 ////////////////////////////////////////////////////////////////////////
107 void AliJetFinder::SetPlotMode(Bool_t b)
108 {
109   // Sets the plotting mode
110   fPlotMode=b;
111   if (b && !fPlots) fPlots = new AliJetControlPlots(); 
112 }
113
114 ////////////////////////////////////////////////////////////////////////
115 TTree* AliJetFinder::MakeTreeJ(char* name)
116 {
117     // Create the tree for reconstructed jets
118     fTreeJ = new TTree(name, "AliJet");
119     fTreeJ->Branch("FoundJet",   &fJets,   1000);
120     fTreeJ->Branch("GenJet",     &fGenJets,1000);
121     fTreeJ->Branch("LeadingPart",&fLeading,1000);
122     return fTreeJ;
123 }
124
125 ////////////////////////////////////////////////////////////////////////
126 void AliJetFinder::WriteRHeaderToFile()
127 {
128   // write reader header
129     AliJetReaderHeader *rh = fReader->GetReaderHeader();
130     rh->Write();
131 }
132
133 ////////////////////////////////////////////////////////////////////////
134 void AliJetFinder::Run()
135 {
136   // Do some initialization
137   Init();
138   // connect files
139   fReader->OpenInputFiles();
140
141   // write headers
142   WriteHeaders();
143   // loop over events
144   Int_t nFirst, nLast, option, debug, arrayInitialised; 
145   nFirst = fReader->GetReaderHeader()->GetFirstEvent();
146   nLast  = fReader->GetReaderHeader()->GetLastEvent();
147
148   option = fReader->GetReaderHeader()->GetDetector();
149   debug  = fReader->GetReaderHeader()->GetDebug();
150   arrayInitialised = fReader->GetArrayInitialised();
151
152   // loop over events
153   for (Int_t i=nFirst;i<nLast;i++) {
154       fReader->FillMomentumArray();
155       fLeading->FindLeading(fReader);
156       fReader->GetGenJets(fGenJets);
157
158       if (option == 0) { // TPC with fMomentumArray
159           if(debug > 1) 
160               printf("In FindJetsC() routine: find jets with fMomentumArray !!!\n");
161           FindJetsC();
162       } else {
163         if(debug > 1) printf("In FindJets() routine: find jets with fUnitArray !!!\n");
164         FindJets();
165       }
166       if (fOut) {
167           fOut->cd();
168       }
169       
170       if (fPlots) fPlots->FillHistos(fJets);
171       fLeading->Reset();
172       fGenJets->ClearJets();
173       Reset();
174   } 
175
176   // write out
177   if (fPlots) {
178       fPlots->Normalize();
179       fPlots->PlotHistos();
180   }
181   if (fOut) {
182       fOut->cd();
183       fPlots->Write();
184       fOut->Close();
185   }
186 }
187
188
189 //
190 // The following methods have been added to allow for event steering from the outside
191 //
192
193 ////////////////////////////////////////////////////////////////////////
194 void AliJetFinder::ConnectTree(TTree* tree, TObject* data)
195 {
196     // Connect the input file
197     fReader->ConnectTree(tree, data);
198 }
199
200 ////////////////////////////////////////////////////////////////////////
201 void AliJetFinder::WriteHeaders()
202 {
203     // Write the Headers
204     TFile* f = new TFile("jets_local.root", "recreate");
205     WriteRHeaderToFile();
206     WriteJHeaderToFile();
207     f->Close();
208 }
209
210 ////////////////////////////////////////////////////////////////////////
211 Bool_t AliJetFinder::ProcessEvent()
212 {
213   //
214   // Process one event
215   // Charged only jets
216   //
217
218   Bool_t ok = fReader->FillMomentumArray();
219   if (!ok) return kFALSE;
220   
221   // Leading particles
222   fLeading->FindLeading(fReader);
223   // Jets
224   FindJets(); // V1
225   //  FindJetsC(); // V2
226   
227   if (fPlots) fPlots->FillHistos(fJets);
228   fLeading->Reset();
229   fGenJets->ClearJets();
230   Reset();  
231   return kTRUE;
232 }
233
234 ////////////////////////////////////////////////////////////////////////
235 Bool_t AliJetFinder::ProcessEvent2()
236 {
237   //
238   // Process one event
239   // Charged only or charged+neutral jets
240   //
241
242   TRefArray* ref = new TRefArray();
243   Bool_t procid = kFALSE;
244   Bool_t ok = fReader->ExecTasks(procid,ref);
245
246   // Delete reference pointer  
247   if (!ok) {delete ref; return kFALSE;}
248   
249   // Leading particles
250   fLeading->FindLeading(fReader);
251   // Jets
252   FindJets();
253   
254   Int_t         nEntRef    = ref->GetEntries();
255   vector<Float_t> vtmp;
256
257   for(Int_t i=0; i<nEntRef; i++)
258     { 
259       // Reset the UnitArray content which were referenced
260       ((AliJetUnitArray*)ref->At(i))->SetUnitTrackID(0);
261       ((AliJetUnitArray*)ref->At(i))->SetUnitEnergy(0.);
262       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag(kPtSmaller);
263       ((AliJetUnitArray*)ref->At(i))->SetUnitCutFlag2(kPtSmaller);
264       ((AliJetUnitArray*)ref->At(i))->SetUnitPxPyPz(kTRUE,vtmp);
265       ((AliJetUnitArray*)ref->At(i))->SetUnitSignalFlag(kBad);
266       ((AliJetUnitArray*)ref->At(i))->SetUnitDetectorFlag(kTpc);
267       ((AliJetUnitArray*)ref->At(i))->SetUnitFlag(kOutJet);
268
269       // Reset process ID
270       AliJetUnitArray* uA = (AliJetUnitArray*)ref->At(i);
271       uA->ResetBit(kIsReferenced);
272       uA->SetUniqueID(0);     
273     }
274
275   // Delete the reference pointer
276   ref->Delete();
277   delete ref;
278
279   if (fPlots) fPlots->FillHistos(fJets);
280   fLeading->Reset();
281   fGenJets->ClearJets();
282   Reset();
283
284   return kTRUE;
285 }
286
287 ////////////////////////////////////////////////////////////////////////
288 void AliJetFinder::FinishRun()
289 {
290     // Finish a run
291     if (fPlots) {
292         fPlots->Normalize();
293         fPlots->PlotHistos();
294     }
295     
296     if (fOut) {
297          fOut->cd();
298          if (fPlots) {
299              fPlots->Write();
300          }
301          fOut->Close();
302     }
303 }
304
305 ////////////////////////////////////////////////////////////////////////
306 void AliJetFinder::AddJet(AliAODJet p)
307 {
308 // Add new jet to the list
309   new ((*fAODjets)[fNAODjets++]) AliAODJet(p);
310 }
311
312 ////////////////////////////////////////////////////////////////////////
313 void AliJetFinder::ConnectAOD(AliAODEvent* aod)
314 {
315 // Connect to the AOD
316     fAODjets = aod->GetJets();
317 }
318
319 ////////////////////////////////////////////////////////////////////////
320 void AliJetFinder::ConnectAODNonStd(AliAODEvent* aod,const char *bname)
321 {
322
323   fAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(bname));
324   // how is this is reset? Cleared?
325 }
326