]> git.uio.no Git - u/mrichter/AliRoot.git/blob - JETAN/AliUA1JetFinder.cxx
Dixed a crash when no proper files/data to display (Mikolaj)
[u/mrichter/AliRoot.git] / JETAN / AliUA1JetFinder.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 // UA1 Jet finder 
18 // manages the search for jets 
19 // Author: jgcn@mda.cinvestav.mx
20 // (code adapted from EMCAL directory)
21 //---------------------------------------------------------------------
22
23 #include <Riostream.h>
24 #include <TLorentzVector.h>
25 #include <TFile.h>
26 #include <TH2F.h>
27 #include "AliUA1JetFinder.h"
28 #include "AliUA1JetHeader.h"
29 #include "UA1Common.h"
30 #include "AliJetReaderHeader.h"
31 #include "AliJetReader.h"
32 #include "AliJet.h"
33
34
35 ClassImp(AliUA1JetFinder)
36
37 ////////////////////////////////////////////////////////////////////////
38
39 AliUA1JetFinder::AliUA1JetFinder():
40   fHeader(0x0),
41   fLego(0x0)
42 {
43   // Default constructor
44 }
45
46 ////////////////////////////////////////////////////////////////////////
47
48 AliUA1JetFinder::~AliUA1JetFinder()
49
50 {
51   // destructor
52 }
53
54 ////////////////////////////////////////////////////////////////////////
55
56 #ifndef WIN32
57 # define ua1_jet_finder ua1_jet_finder_
58 # define hf1 hf1_
59 # define type_of_call
60
61 #else
62 # define ua1_jet_finder UA1_JET_FINDER
63 # define hf1 HF1
64 # define type_of_call _stdcall
65 #endif
66
67 extern "C" void type_of_call
68 ua1_jet_finder(Int_t& ncell, Int_t& ncell_tot,
69                Float_t  etc[60000],  Float_t etac[60000],
70                Float_t  phic[60000],
71                Float_t& min_move, Float_t& max_move, Int_t& mode,
72                Float_t& prec_bg,  Int_t& ierror);
73
74 extern "C" void type_of_call hf1(Int_t& id, Float_t& x, Float_t& wgt);
75
76 ////////////////////////////////////////////////////////////////////////
77
78 void AliUA1JetFinder::FindJetsTPC()
79
80 {
81   // initialize event, look for jets, download jet info
82
83   // initialization in 2 steps
84   // 1) transform input to pt,eta,phi plus lego
85   // 2) dump lego
86   
87   // 1) transform input to pt,eta,phi plus lego
88   TClonesArray *lvArray = fReader->GetMomentumArray();
89   Int_t nIn =  lvArray->GetEntries();
90   if (nIn == 0) return;
91     
92   // local arrays for input
93   Float_t* enT  = new Float_t[nIn];
94   Float_t* ptT  = new Float_t[nIn];
95   Float_t* etaT = new Float_t[nIn];
96   Float_t* phiT = new Float_t[nIn];
97
98   // load input vectors
99   for (Int_t i = 0; i < nIn; i++){
100     TLorentzVector *lv = (TLorentzVector*) lvArray->At(i);
101     enT[i]  = lv->Energy();
102     ptT[i]  = lv->Pt();
103     etaT[i] = lv->Eta();
104     phiT[i] = ((lv->Phi() < 0) ? (lv->Phi()) + 2 * TMath::Pi() : lv->Phi());
105     if (fReader->GetCutFlag(i) == 1) 
106       fLego->Fill(etaT[i], phiT[i], ptT[i]);
107   }
108   fJets->SetNinput(nIn);
109   
110   // 2) dump lego
111   // check enough space! *to be done*
112   Float_t etCell[60000];   //! Cell Energy
113   Float_t etaCell[60000];  //! Cell eta
114   Float_t phiCell[60000];  //! Cell phi
115
116   Int_t nCell = 0;
117   TAxis* xaxis = fLego->GetXaxis();
118   TAxis* yaxis = fLego->GetYaxis();
119   Float_t e = 0.0;
120   for (Int_t i = 1; i <= fHeader->GetLegoNbinEta(); i++) {
121       for (Int_t j = 1; j <= fHeader->GetLegoNbinPhi(); j++) {
122           e = fLego->GetBinContent(i,j);
123           if (e > 0.0) e -= fHeader->GetMinCellEt();
124           if (e < 0.0) e = 0.;
125           Float_t eta  = xaxis->GetBinCenter(i);
126           Float_t phi  = yaxis->GetBinCenter(j);            
127           etCell[nCell]  = e;
128           etaCell[nCell] = eta;
129           phiCell[nCell] = phi;
130           nCell++;
131       }
132   }
133
134   // run the algo. Parameters from header
135   Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
136   Float_t minmove = fHeader->GetMinMove();
137   Float_t maxmove = fHeader->GetMaxMove();
138   Int_t mode      = fHeader->GetMode();
139   Float_t precbg  = fHeader->GetPrecBg();
140   Int_t ierror;
141   ua1_jet_finder(nCell, nTot, etCell, etaCell, phiCell, 
142                  minmove, maxmove, mode, precbg, ierror);
143
144   // sort jets
145   Int_t * idx  = new Int_t[UA1JETS.njet];
146   TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
147
148   // download jet info.   
149   for(Int_t i = 0; i < UA1JETS.njet; i++) {
150     // reject events outside acceptable eta range
151     if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
152         || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
153           continue;
154       Float_t px, py,pz,en; // convert to 4-vector
155       px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
156       py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
157       pz = UA1JETS.etj[idx[i]] /
158           TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
159       en = TMath::Sqrt(px * px + py * py + pz * pz);
160       fJets->AddJet(px, py, pz, en);
161   }
162   
163   // find multiplicities and relationship jet-particle
164   // find also percentage of pt from pythia
165   Int_t* injet = new Int_t[nIn];
166   Int_t* sflag = new Int_t[nIn];
167   for (Int_t i = 0; i < nIn; i++) {injet[i]= 0;sflag[i]=0;}
168   Int_t* mult  = new Int_t[fJets->GetNJets()];
169   Int_t* ncell  = new Int_t[fJets->GetNJets()];
170   Float_t* percentage  = new Float_t[fJets->GetNJets()];
171
172   for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
173       Float_t pt_sig = 0.0;
174       mult[i] = 0;
175       ncell[i] = UA1JETS.ncellj[i];
176       for (Int_t j = 0; j < nIn; j++) {
177           Float_t deta = etaT[j] - fJets->GetEta(i);
178           Float_t dphi = phiT[j] - fJets->GetPhi(i);
179           if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
180           if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
181           Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
182           if (dr < fHeader->GetRadius() && injet[j] == 0) {
183               injet[j] = -(i+1);
184               if ((fReader->GetCutFlag(j)) == 1 &&
185                   (etaT[j] < fHeader->GetLegoEtaMax()) &&
186                   (etaT[j] > fHeader->GetLegoEtaMin())) {
187                   injet[j] = i+1;
188                   mult[i]++;
189                   if (fReader->GetSignalFlag(j) == 1) {
190                     pt_sig+=ptT[j];
191                     sflag[j]=1;
192                   }
193               }
194           }
195       }
196       percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
197         ((Double_t) fJets->GetPt(i));    
198   }
199   fJets->SetNCells(ncell);
200   fJets->SetPtFromSignal(percentage);
201   fJets->SetMultiplicities(mult);
202   fJets->SetInJet(injet);
203   fJets->SetEtaIn(etaT);
204   fJets->SetPhiIn(phiT);
205   fJets->SetPtIn(ptT);
206   fJets->SetEtAvg(UA1JETS.etavg);
207   delete ncell;
208   delete enT;
209   delete ptT;
210   delete etaT;
211   delete phiT;
212   delete injet;
213   delete idx;
214   delete mult;
215   delete percentage;
216 }
217
218 void AliUA1JetFinder::FindJets()
219 { // Find jets with TPC or EMCAL or TPC + EMCAL information
220   // initialize event, look for jets, download jet info
221
222   // 1) transform input to pt,eta,phi plus lego
223     AliJetUnitArray *fUnit = fReader->GetUnitArray();
224     Int_t nIn = fUnit->GetUnitEntries();
225     Int_t fOption = fReader->GetReaderHeader()->GetDetector();
226     Int_t fDebug = fReader->GetReaderHeader()->GetDebug();
227     
228     if(fDebug>1){ 
229         printf("Inside FindJets function ! \n");
230         printf("fOption : %d \n",fOption);
231     }
232     
233     if(fDebug>10){ 
234         cout << "fUnit : " << fUnit << endl;
235         printf("nIn = fUnit->GetUnitEntries() : %d \n",nIn);
236     }
237     
238     if(fDebug>10){
239         printf("     -----------------------------------------------------------------\n");
240         printf("     --- All inputs in fUnitArray ---\n");
241         for(Int_t i=0; i<nIn ; i++){
242             if(fUnit[i].GetUnitEnergy()!=0){
243                 printf("     -----------------------------------------------------------------\n");
244                 cout << "     |  ID : " << fUnit[i].GetUnitID() << "  |  Eta : " << fUnit[i].GetUnitEta() << "  |  Phi : " << fUnit[i].GetUnitPhi() << "  |  Energy : " << fUnit[i].GetUnitEnergy() << " |" <<endl;
245             }
246         }
247     }
248     
249     if (nIn == 0) return;
250     
251     Int_t nCandidateCut = 0;
252     Int_t nCandidateCut2 = 0;
253     Int_t nCandidate = 0;
254     for (Int_t i = 0; i < nIn; i++){
255         if(fUnit[i].GetUnitEnergy()>0. && fUnit[i].GetUnitCutFlag()==1) {
256             // Number of good tracks/digits which have passed the pt cut
257             nCandidateCut += 1;
258         }
259     if(fUnit[i].GetUnitEnergy()>0.) {
260         // Number of good tracks/digits without pt cut
261         nCandidate += 1;
262     }
263     }
264     
265     if(fDebug>=3){
266     cout << "nCandidate : " << nCandidate << endl;
267     cout << "nCandidateCut : " << nCandidateCut << endl;
268     cout << "nCandidateCut2 : " << nCandidateCut2 << endl;
269     }
270     
271   // local arrays for input No Cuts
272   // Both pt < ptMin and pt > ptMin
273     Float_t* enT       = new Float_t[nCandidate];
274     Float_t* ptT       = new Float_t[nCandidate];
275     Float_t* etaT      = new Float_t[nCandidate];
276     Float_t* phiT      = new Float_t[nCandidate];
277     Float_t* detaT     = new Float_t[nCandidate];
278     Float_t* dphiT     = new Float_t[nCandidate];
279     Float_t* cFlagT    = new Float_t[nCandidate];
280     Float_t* cClusterT = new Float_t[nCandidate];
281     Float_t* idT       = new Float_t[nCandidate];
282     Int_t loop1 = 0;
283     
284     fJets->SetNinput(nCandidate);
285     
286     if(fDebug>3){
287         cout << "nCandidate : " << nCandidate << endl;
288         //    cout << "nMultCandidate : " << nMultCandidate << endl;
289     }
290     
291     Float_t *etCell = new Float_t[nIn];   //! Cell Energy - Extracted from UnitArray
292     Float_t *etaCell = new Float_t[nIn];  //! Cell eta - Extracted from UnitArray
293     Float_t *phiCell = new Float_t[nIn];  //! Cell phi - Extracted from UnitArray
294     
295     // Information extracted from fUnitArray
296     for(Int_t i=0; i<nIn; i++) 
297     {
298         if(fUnit[i].GetUnitCutFlag()==1){ 
299             etCell[i] = fUnit[i].GetUnitEnergy();
300             if (etCell[i] > 0.0) etCell[i] -= fHeader->GetMinCellEt();
301             if (etCell[i] < 0.0) etCell[i] = 0.;
302             etaCell[i] = fUnit[i].GetUnitEta();
303             phiCell[i] = fUnit[i].GetUnitPhi();
304         }
305         else {
306             etCell[i]  = 0.;
307             etaCell[i] = fUnit[i].GetUnitEta();
308             phiCell[i] = fUnit[i].GetUnitPhi();
309         }
310         
311         if(fUnit[i].GetUnitEnergy()>0.){
312             ptT[loop1]   = fUnit[i].GetUnitEnergy();
313             enT[loop1]   = fUnit[i].GetUnitEnergy();
314             etaT[loop1]  = fUnit[i].GetUnitEta();
315             phiT[loop1]  = fUnit[i].GetUnitPhi();
316             detaT[loop1] = fUnit[i].GetUnitDeta();
317             dphiT[loop1] = fUnit[i].GetUnitDphi();
318             cFlagT[loop1]= fUnit[i].GetUnitCutFlag();
319             idT[loop1]   = fUnit[i].GetUnitID();
320             loop1++;
321         }
322     }
323     
324     
325     if(fDebug > 40) // For comparison
326     {
327         for(Int_t j=0; j<nIn; j++) {
328             if(etCell[j]>0){
329                 cout << "etCell[" << j << "] : " << etCell[j] << endl;
330                 cout << "etaCell[" << j << "] : " << etaCell[j] << endl;
331                 cout << "phiCell[" << j << "] : " << phiCell[j] << endl;
332             }
333         }
334     }
335     
336     // Run the algo. Parameters from header
337     //  Int_t nTot      = (fHeader->GetLegoNbinEta())*(fHeader->GetLegoNbinPhi());
338     Int_t nTot      = nIn;
339     Float_t minmove = fHeader->GetMinMove();
340     Float_t maxmove = fHeader->GetMaxMove();
341     Int_t mode      = fHeader->GetMode();
342     Float_t precbg  = fHeader->GetPrecBg();
343     Int_t ierror;
344     
345     ua1_jet_finder(nIn, nTot, etCell, etaCell, phiCell, 
346                    minmove, maxmove, mode, precbg, ierror);
347     
348     // sort jets
349     Int_t * idx  = new Int_t[UA1JETS.njet];
350     TMath::Sort(UA1JETS.njet, UA1JETS.etj, idx);
351     
352     if(fDebug > 20)
353     {
354         for(Int_t i = 0; i < UA1JETS.njet; i++) 
355         {
356             cout << "Number of jets found, UA1JETS.njet : " << UA1JETS.njet << endl;
357             cout << "UA1JETS.etj : " << UA1JETS.etj << endl;
358             cout << "idx[" << i << "] : " << idx[i] << endl;
359             cout << "UA1JETS.etaj[1][" << idx[i] << "] : " << UA1JETS.etaj[1][idx[i]] << endl;
360             cout << "UA1JETS.phij[1][" << idx[i] << "] : " << UA1JETS.phij[1][idx[i]] << endl;
361             cout << "UA1JETS.etj[" << idx[i] << "] : " << UA1JETS.etj[idx[i]] << endl;
362         }
363     }
364     
365     // download jet info.   
366     for(Int_t i = 0; i < UA1JETS.njet; i++) {
367         // reject events outside acceptable eta range
368         if (((UA1JETS.etaj[1][idx[i]])> (fHeader->GetJetEtaMax()))
369             || ((UA1JETS.etaj[1][idx[i]]) < (fHeader->GetJetEtaMin())))
370         {
371             continue;
372         }
373         
374         Float_t px, py,pz,en,pT; // convert to 4-vector
375         px = UA1JETS.etj[idx[i]] * TMath::Cos(UA1JETS.phij[1][idx[i]]);
376         py = UA1JETS.etj[idx[i]] * TMath::Sin(UA1JETS.phij[1][idx[i]]);
377         pz = UA1JETS.etj[idx[i]] /
378             TMath::Tan(2.0 * TMath::ATan(TMath::Exp(-UA1JETS.etaj[1][idx[i]])));
379         en = TMath::Sqrt(px * px + py * py + pz * pz);
380         pT = TMath::Sqrt(px * px + py * py);
381         
382         fJets->AddJet(px, py, pz, en);
383         
384     }
385
386     // find multiplicities and relationship jet-particle
387     // find also percentage of pt from pythia
388     
389     Int_t* injet = new Int_t[nCandidate];
390     Int_t* sflag = new Int_t[nCandidate];
391     for (Int_t i = 0; i < nCandidate; i++) {injet[i]= 0;sflag[i]=0;}
392     Int_t* mult  = new Int_t[fJets->GetNJets()];
393     Int_t* ncell  = new Int_t[fJets->GetNJets()];
394     Float_t* percentage  = new Float_t[fJets->GetNJets()];
395     
396     // Instead of using etaT below, it would be interesting to use the previous fUnitArray object
397     // With the particle ID, it is possible to to have access to its physical properties and one can,
398     // for example, set if the corresponding particle is inside or outside the jet with the flag 
399     // kOutJet/kInJet, other possibilities...
400     
401     for(Int_t i = 0; i < (fJets->GetNJets()); i++) {
402         Float_t pt_sig = 0.0;
403         mult[i] = 0;
404         ncell[i] = UA1JETS.ncellj[i];
405         for (Int_t j = 0; j < nCandidate; j++) {
406             Float_t deta = etaT[j] - fJets->GetEta(i);
407             Float_t dphi = phiT[j] - fJets->GetPhi(i);
408             if (dphi < -TMath::Pi()) dphi= -dphi - 2.0 * TMath::Pi();
409             if (dphi > TMath::Pi()) dphi = 2.0 * TMath::Pi() - dphi;
410             Float_t dr = TMath::Sqrt(deta * deta + dphi * dphi);
411             if (dr < fHeader->GetRadius() && injet[j] == 0) {
412                 injet[j] = -(i+1);
413                 if(cFlagT[j] == 1 &&
414                    (etaT[j] < fHeader->GetLegoEtaMax()) &&
415                    (etaT[j] > fHeader->GetLegoEtaMin())) {
416                     injet[j] = i+1;
417                     mult[i]++;
418                     pt_sig+=enT[j];
419                     sflag[j]=1;
420                 }
421             }
422             if(fDebug>10){
423                 cout << "mult[" << i << "] : " << mult[i] << endl;
424                 cout << "ncell[" << i << "] : " << ncell[i] << endl;
425             }
426         }
427         percentage[i] = (pt_sig-ncell[i]*UA1JETS.etavg)/
428             ((Double_t) fJets->GetPt(i));    
429     }
430     
431     fJets->SetNCells(ncell);
432     fJets->SetPtFromSignal(percentage);
433     fJets->SetMultiplicities(mult);
434     fJets->SetInJet(injet);
435     fJets->SetEtaIn(etaT);
436     if(fDebug>10){
437         for(Int_t i=0; i<nCandidate ; i++){
438             cout << "phiT[" << i << "] : " << phiT[i] << endl;
439             cout << "etaT[" << i << "] : " << etaT[i] << endl;
440         }
441     }
442     fJets->SetPhiIn(phiT);
443     fJets->SetPtIn(enT);
444     fJets->SetEtAvg(UA1JETS.etavg);
445     delete etCell;
446     delete etaCell;
447     delete phiCell;
448     delete ncell;
449     delete cFlagT;
450     delete cClusterT;
451     delete enT;
452     delete ptT;
453     delete etaT;
454     delete phiT;
455     delete detaT;
456     delete dphiT;
457     delete injet;
458     delete idx;
459     delete mult;
460     delete percentage;
461     
462 }
463
464
465 ////////////////////////////////////////////////////////////////////////
466
467 void AliUA1JetFinder::Reset()
468 {
469   fLego->Reset();
470   fJets->ClearJets();
471 }
472
473 ////////////////////////////////////////////////////////////////////////
474
475 void AliUA1JetFinder::WriteJHeaderToFile()
476 {
477   fOut->cd();
478   fHeader->Write();
479 }
480
481 ////////////////////////////////////////////////////////////////////////
482
483 void AliUA1JetFinder::Init()
484 {
485   // initializes some variables
486   Float_t dEta, dPhi;
487   dEta = (fHeader->GetLegoEtaMax()-fHeader->GetLegoEtaMin())
488     /((Float_t) fHeader->GetLegoNbinEta());
489   dPhi = (fHeader->GetLegoPhiMax()-fHeader->GetLegoPhiMin())
490     /((Float_t) fHeader->GetLegoNbinPhi());
491
492   UA1CELL.etaCellSize = dEta;
493   UA1CELL.phiCellSize = dPhi;
494
495   // jet parameters
496   UA1PARA.coneRad = fHeader->GetRadius();
497   UA1PARA.etSeed  = fHeader->GetEtSeed();
498   UA1PARA.ejMin   = fHeader->GetMinJetEt();
499   UA1PARA.etMin   = fHeader->GetMinCellEt();
500
501   // book lego
502   fLego = new 
503     TH2F("legoH","eta-phi",
504          fHeader->GetLegoNbinEta(), fHeader->GetLegoEtaMin(), 
505          fHeader->GetLegoEtaMax(),  fHeader->GetLegoNbinPhi(), 
506          fHeader->GetLegoPhiMin(), fHeader->GetLegoPhiMax());
507 }