1 //*-- Author: Aleksei Pavlinov (WSU)
4 #if !defined(__CINT__) || defined(__MAKECINT__)
5 // These files are needed for compilation only
6 // (see directory /auto/u/pavlinov/macros on PDSF)
7 // Using with caution if you are not PAI
8 #include "macroIncludePai.h"
9 #include "macroIncludeAlice.h"
10 extern "C++" {void loadlibs();}
13 // For convinience only
14 class AliEMCALJetFinder; AliEMCALJetFinder* jetFinder=0;
15 TH2F *hLego, *hLegoW, *hTrig;
16 TH1F *hETrigPatch, *hJetPartPt;
18 // directory with file and number of files
19 TString dirIn("/auto/alice/sahal/PPR3/11/"), nFileIn;
20 Int_t file1=1, file2=20, nevMax=0;
21 // change dirOut if you want.
22 TString dirOut("RES/FILE/"), nFileOut;
24 // Default value for jet fineer - 11-feb-2202
26 Float_t coneRadius = 0.3;
28 Float_t minJetEt = 40.;
29 Float_t minCellEt = 1.; // be carefull
30 Float_t minPtCut = 2.; // for including charge particles to file
32 Int_t modeSubtraction = 1;
33 Float_t minMove = 0.05;
34 Float_t maxMove = 0.15;
35 Float_t precBg = 0.035;
37 Bool_t hadronCorrection = kTRUE;
39 class AliEMCALJetMicroDst; AliEMCALJetMicroDst *dst=0;
40 Int_t nevSum = 0, nevF;
41 TArrayI nevInFiles(100);
44 void jetDst(Int_t mode, Int_t var, Int_t nmax)
47 // mode - variable for selection input mode (see function defineInput());
48 // var - variable for selection jetfinder mode (see function defineJetFinder()).
50 // Dynamically link some shared libs
51 if (gClassTable->GetID("AliRun") < 0) {
52 gROOT->LoadMacro("../macros/loadlibs.C");
56 if(!defineInput(mode)) return;
58 dst = new AliEMCALJetMicroDst;
60 if(!nFileOut.Contains(".root")){
64 dst->Create(nFileOut.Data());
65 gROOT->GetListOfBrowsables()->Add(dst);
67 for(Int_t j=file1; j<=file2; j++){
69 if(file2>1) {// multiple file
71 if (mode>=11 && mode<=11) {
73 nFileIn += ("/galice.root");
74 } else if(mode>=21 && mode<=51){
80 // Connect the Root Galice file containing Geometry, Kine and Hits
81 file = new TFile(nFileIn.Data(), "READ");
82 printf("\n FILE %i -> %s\n", j, nFileIn.Data());
84 if (file == 0) continue;
86 // Get AliRun object from file or create it if not on file
88 gAlice = (AliRun*)(file->Get("gAlice"));
90 printf("AliRun object found on file\n");
92 defineJetFinder(mode,var);
94 // Loop over events in file
97 // if(!c1) c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
100 if(nmax>0) nevMax = nmax;
101 for (Int_t nev = 0; nev<= nevMax; nev++) {
102 jetFinder->ResetJets();
103 Int_t nparticles = gAlice->GetEvent(nev);
104 if (nparticles <= 0) break; // last event
109 jetFinder->FillFromPartons();
111 for(Int_t ij=0; ij<jetFinder->Njets(); ij++){
112 hJetPartPt->Fill(jetFinder->JetEnergy(ij));
114 jetFinder->ResetJets();
117 jetFinder->FillFromHits();
119 // gROOT->cd(); be carefull with this staff
120 hLegoW = new TH2F((*jetFinder->GetLego())); // for saving
121 sw = hLegoW->GetName();
123 hLegoW->SetName(sw.Data());
124 sw = "Energy in ECAL:Event";
126 hLegoW->SetTitle(sw.Data());
129 hLego = jetFinder->GetLego();
130 hTrig = bookTrigHist(hLego);
132 fillTriggerPatch(hLegoW, hTrig, hETrigPatch);
134 // jetFinder->FillFromDigits();
136 jetFinder->FillFromTracks(1, 0);
137 // ^ preserves info from hit
139 // TPC information from Hits associated to the EMCAL
140 // jetFinder->FillFromHitFlaggedTracks(1);
143 dst->Fill(gAlice, jetFinder);
146 printf("\n Events : %d Jets : %d \n", nev, jetFinder->Njets());
147 Int_t njet = jetFinder->Njets();
148 for (Int_t nj=0; nj<njet; nj++)
150 printf("\n Jet Energy %5d %8.2f \n",
151 nj, jetFinder->JetEnergy(nj));
152 printf("Jet Phi %5d %8.2f %8.2f \n",
153 nj, jetFinder->JetPhiL(nj), jetFinder->JetPhiW(nj));
154 printf("Jet Eta %5d %8.2f %8.2f \n",
155 nj, jetFinder->JetEtaL(nj), jetFinder->JetEtaW(nj));
159 if (c1 && nev <-10 && njet==1 && jetFinder->JetEnergy(0)>150) {
161 cout << "Event "<< nev <<endl;
168 nevInFiles[j-1] = nevF;
170 delete gAlice; // don't delete after closing file
174 if(c1) hLego->Draw("lego");
176 if(nmax==0) dst->Close();
177 nevInFiles.Set(file2);
178 printf("\n Summ.events %6i -> %6i\n", nevSum, (Int_t) nevInFiles.GetSum());
181 for(Int_t j=file1; j<=file2; j++){
182 printf(" %i4 ",nevInFiles[j-1]);
183 if((j-file1)%10==0 && j!=file1) printf("\n ");
187 // ###################################################################
188 TString* defineInput(Int_t mode)
192 // case -5 needing for testing after merging
195 hadronCorrection = kFALSE; // for checking Tom's idea
199 if (mode==-2) minPtCut = 1.;
200 else if(mode<=-3 || mode>=-5) {
201 minCellEt = 0.0; // 13-feb-2002
213 nFileIn = "galice.root";
214 sprintf(ctmp,"PythiaMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
215 mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
216 nFileOut = dirOut + ctmp;
217 if(!hadronCorrection) nFileOut.ReplaceAll(".root","_NoHadrCorr.root");
220 nFileOut = dirOut + "dstTest.root";
224 // jj100 41 events per file - see Log file
225 dirIn = "/auto/alice/sahal/PPR4/1/";
227 file2 = 30; // max 75
229 sprintf(ctmp,"HijingMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
230 mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
231 nFileOut = dirOut + ctmp;
233 //=============================================================
234 // central events with jet trigger - 15 feb 2002
235 //=============================================================
237 dirIn = "/auto/alice/sahal/kHijing_jj50/";
241 sprintf(ctmp,"Hijing_jj50Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
242 mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
243 nFileOut = dirOut + ctmp;
246 dirIn = "/auto/alice/sahal/kHijing_jj75/";
250 sprintf(ctmp,"Hijing_jj75Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
251 mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
252 nFileOut = dirOut + ctmp;
255 dirIn = "/auto/alice/sahal/kHijing_jj100/";
259 sprintf(ctmp,"Hijing_jj100Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
260 mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
261 nFileOut = dirOut + ctmp;
264 dirIn = "/auto/alice/sahal/kHijing_jj125/";
268 sprintf(ctmp,"Hijing_jj125Mode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
269 mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
270 nFileOut = dirOut + ctmp;
274 // BG - HIJING events
275 dirIn = "/auto/alice/sahal/background/";
279 sprintf(ctmp,"HijingBGMode%iR%3.1fSeed%3.1fJEt%4.1fCEt%3.1fPtCut%3.1f.root",
280 mode, coneRadius, etSeed, minJetEt, minCellEt, minPtCut);
281 nFileOut = dirOut + ctmp;
284 printf("\n<I> Wrong input mode => %i\ncout", mode);
287 printf("<I> dir %s file1 %3i file2 %3i\n", dirIn.Data(), file1, file2);
291 void defineJetFinder(Int_t mode, Int_t var)
293 // Create and configure jet finder
294 // gAlice must be defined already !!!
296 jetFinder = new AliEMCALJetFinder("UA1 Jet Finder", "Pre Production");
299 jetFinder->SetDebug(debugJet);
300 jetFinder->SetConeRadius(coneRadius);
301 jetFinder->SetEtSeed(etSeed);
302 jetFinder->SetMinJetEt(minJetEt);
303 jetFinder->SetMinCellEt(minCellEt);
304 jetFinder->SetPtCut(minPtCut);
306 if(var>=2) jetFinder->SetMomentumSmearing(1);
307 if(var>=3) jetFinder->SetEfficiencySim(1);
309 jetFinder->SetHadronCorrection(0);
310 if(var>=4 && hadronCorrection) {
311 jetFinder->SetHadronCorrection(1);
312 jetFinder->SetHadronCorrector(AliEMCALHadronCorrectionv0::Instance());
314 if(var>=5) { // 11-feb-2002
315 jetFinder->SetParametersForBgSubtraction
316 (modeSubtraction,minMove,maxMove,precBg);
319 if(mode<0) jetFinder->SetSamplingFraction(1.);//must use only for PYTHIA file
320 else jetFinder->SetSamplingFraction(12.9);
321 gROOT->GetListOfBrowsables()->Add(jetFinder);
325 void testDST(char* fname)
327 dst = new AliEMCALJetMicroDst;
330 gROOT->GetListOfBrowsables()->Add(dst);
333 //----------------- 4-feb-2002 --------------
334 TH2F* bookTrigHist(TH2F *hid)
337 TAxis *xax = hid->GetXaxis();
338 Double_t xmin = xax->GetXmin();
339 Double_t xmax = xax->GetXmax();
341 TAxis *yax = hid->GetYaxis();
342 Double_t ymin = yax->GetXmin();
343 Double_t ymax = yax->GetXmax();
345 TH2F *hout = new TH2F("hTrig","Energy grid for trigger patch",
346 nx,xmin,xmax, ny,ymin,ymax);
347 hout->AddDirectory(0);
350 hETrigPatch = new TH1F("hETrigPatch","Energy in trigger patch",300, 0.0, 300.);
351 hJetPartPt = new TH1F("hJetPartPt", "E_{T} for partons jet", 300, 0.0, 300.);
352 hETrigPatch->AddDirectory(0);
356 void fillTriggerPatch(TH2F *hid, TH2F *htrig, TH1F *hE)
358 TAxis *xax = hid->GetXaxis();
359 TAxis *yax = hid->GetYaxis();
360 Int_t nx = xax->GetNbins();
361 Int_t ny = yax->GetNbins();
365 for(Int_t ix=1; ix<=nx; ix++){
366 x = xax->GetBinCenter(ix);
367 for(Int_t iy=1; iy<=ny; iy++){
368 y = yax->GetBinCenter(iy);
369 e = hid->GetBinContent(ix, iy);
374 xax = htrig->GetXaxis();
375 yax = htrig->GetYaxis();
376 nx = xax->GetNbins();
377 ny = yax->GetNbins();
378 for(Int_t ix=1; ix<=nx; ix++){
379 for(Int_t iy=1; iy<=ny; iy++){
380 hE->Fill(htrig->GetBinContent(ix, iy));