]>
Commit | Line | Data |
---|---|---|
498224e8 | 1 | //*-- Author: Aleksei Pavlinov(WSU) |
2 | // 8-feb-2002 new version with AliEMCALJetMicroDst | |
3 | #include "anaDst.h" | |
4 | #include "jetDst.h" | |
5 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
6 | // These files are needed for compilation only | |
7 | // (see directory /auto/u/pavlinov/macros on PDSF) | |
8 | // Using with caution if you are not PAI | |
9 | #include "macroIncludePai.h" | |
10 | #include "macroIncludeAlice.h" | |
30bd7436 | 11 | extern "C++" {void loadlibs();} |
498224e8 | 12 | #endif |
13 | ||
14 | class TH1F; | |
15 | TH1F* hJetEt, *hJetEta, *hJetPhi, *hNtracksInJet, *hPtTracksInJet, *hTrackDR, *hNjet; | |
16 | ||
17 | class TH2F; TH2F* hPartPt, *hPartEta, *hPartPhi; | |
18 | TH1F* hPartDiffPt, *hPartDiffEta, * hPartDiffPhi; | |
19 | ||
20 | TH1F *hJetPartAng, *hDeta, *hDphi, *hDpt; // Comparing between jet and nearest parton | |
21 | ||
22 | class TCanvas; TCanvas *c1, *c2, *c3; | |
23 | class TPad; TPad *p1; | |
24 | class TList; TList *lPartons, *lJets, *lCompPJ; | |
25 | ||
26 | class TVector3; | |
27 | ||
28 | Int_t nparticles, nJet; | |
29 | Double_t phiMax = TMath::Pi()*120./180.; | |
30 | Double_t rJet = 0.5; // this default | |
31 | ||
32 | TFile *file; | |
33 | class AliEMCALJetMicroDst; AliEMCALJetMicroDst *jDst=0; | |
34 | Int_t nevMax = 1000; | |
35 | TString fname; | |
36 | // fitting | |
37 | class TF1 *g1;// gauss | |
498224e8 | 38 | |
39 | void | |
40 | anaDst(Int_t mode) | |
41 | { | |
42 | // Dynamically link some shared libs | |
43 | if (gClassTable->GetID("AliRun") < 0) { | |
44 | gROOT->LoadMacro("../macros/loadlibs.C"); | |
45 | loadlibs(); | |
46 | } | |
47 | ||
48 | //partons after HSC and jet in acceptance. | |
49 | TVector3 *vJet; | |
50 | TVector3 *vPart1, *vPart2; | |
51 | vJet = new TVector3; | |
52 | vPart1 = new TVector3; | |
53 | vPart2 = new TVector3; | |
54 | ||
55 | // for getting information about nFileOut | |
56 | //#if !defined(__CINT__) | |
57 | // gROOT->ProcessLine(".L jetDst.C+"); | |
58 | //#else | |
59 | // gROOT->LoadMacro("jetDst.C"); | |
60 | //#endif | |
61 | if(!defineMicroDst(mode)) return; | |
62 | ||
63 | bookHist1(); | |
64 | bookHistPartonsAfterHardSc(); | |
65 | ||
66 | Float_t phiT[50], etaT[50], ptT[50]; | |
67 | ||
68 | Int_t nbytes=0, nb=0, nentries = Int_t(jDst->GetTree()->GetEntries()); | |
69 | if(nevMax>nentries) nevMax = nentries; | |
70 | for (int nev=0; nev< nevMax; nev++) { | |
71 | nb = jDst->GetEntry(nev); | |
72 | nbytes += nb; | |
73 | printf("\n Event .............%d bytes %d\n", nev, nb); | |
74 | // jDst->Print(); | |
75 | if(nb<0) break; // last event or something wrong | |
76 | // return; | |
77 | ||
78 | Bool_t inEmc = fillInfoForPartons(vPart1, vPart2); | |
79 | // continue; | |
80 | ||
81 | nJet = jDst->GetNjet(); | |
82 | Double_t et, eta, phi; | |
83 | ||
84 | for(Int_t ij=0; ij < nJet; ij++) { | |
85 | ||
86 | jDst->GetJet(ij, 1, (*vJet)); // W | |
87 | ||
88 | et = vJet->Pt(); | |
89 | eta = vJet->Eta(); | |
90 | phi = vJet->Phi(); | |
91 | ||
92 | printf("\n Jet:%d E %f phi %f eta %f tracks %d\n", ij, | |
93 | et, phi, eta, 0); | |
94 | ||
95 | // hNtracksInJet->Fill((Float_t)lJet->NTracks()); | |
96 | // Safe zone - rJet - 13-feb-2002 | |
97 | if(nJet==1 && TMath::Abs(eta)<(0.7-rJet)&& | |
98 | phi>rJet && phi<(phiMax-rJet)) { | |
99 | hJetEt->Fill(et); | |
100 | hJetEta->Fill(eta); | |
101 | hJetPhi->Fill(phi); | |
102 | if(inEmc) compareJetPartons(vJet, vPart1, vPart2); | |
103 | } | |
104 | ||
105 | } // jet | |
106 | cleanUpEvent(); | |
107 | } // event | |
108 | ||
109 | gStyle->SetOptStat(1111111); | |
110 | ||
111 | //drawPartons(); | |
112 | ||
113 | #if !defined(__CINT__) | |
114 | if(hJetPartAng&&hJetPartAng->GetEntries()>2.){ | |
115 | _drawListOfHist(lCompPJ, 2, 2); | |
116 | fitGauss(hDpt, 0); | |
117 | } | |
118 | #endif | |
119 | ||
120 | gROOT->cd(); | |
121 | // printf("\n Variant %1i : Alice Root file => %s \n", | |
122 | //mode, nf.Data()); | |
123 | } | |
124 | ||
125 | Bool_t defineMicroDst(Int_t mode) | |
126 | { | |
127 | TString dir("/auto/u/pavlinov/ALICE/RES/FILE/"); | |
128 | // TString *sTmp = defineInput(mode); // see jetDst.h | |
129 | switch(mode){ | |
130 | // case -4: | |
131 | //if(sTmp) fname = sTmp->Data(); | |
132 | //else return kFALSE; | |
133 | //break; | |
134 | case -3: | |
135 | fname = dir + "PythiaMode-3R0.5Seed4.0JEt20.0CEt0.0PtCut0.0.root"; | |
136 | break; | |
137 | case -2: | |
138 | fname = dir + "PythiaMode-2R0.5Seed4.0JEt20.0CEt1.0PtCut1.0.root"; | |
139 | break; | |
140 | case -1: | |
141 | fname = dir + "PythiaMode-1R0.5Seed4.0JEt20.0CEt1.0PtCut2.0.root"; | |
142 | break; | |
143 | ||
144 | case 11: | |
145 | fname = dir + "HijingMode11R0.3Seed4.0JEt40.0CEt1.0PtCut2.0.root"; | |
146 | break; | |
147 | default: | |
148 | printf("\n<I> Wrong input mode => %i\n", mode); | |
149 | return kFALSE; | |
150 | } | |
151 | jDst = new AliEMCALJetMicroDst; | |
152 | jDst->Open(fname.Data()); | |
153 | gROOT->GetListOfBrowsables()->Add(jDst); | |
154 | return kTRUE; | |
155 | } | |
156 | ||
157 | Bool_t fillInfoForPartons(TVector3 *vec1, TVector3 *vec2) | |
158 | { | |
159 | // npart=4 for HIJING and 2 for PYTHIA | |
160 | Int_t npart = jDst->GetNpart(); | |
161 | if(!jDst->GetParton(npart-2, (*vec1))||!jDst->GetParton(npart-1, (*vec2)) ) | |
162 | return kFALSE; | |
163 | ||
164 | static Double_t dphi; | |
165 | ||
166 | Bool_t inEmc = kFALSE; | |
167 | TVector3 *vec = vec1; | |
168 | for (Int_t i=0; i<2; i++) { | |
169 | Double_t ri = Double_t(i); | |
170 | hPartPt->Fill(vec->Pt(), ri); | |
171 | hPartEta->Fill(vec->Eta(), ri); | |
172 | hPartPhi->Fill(vec->Phi(), ri); | |
173 | if(fabs(vec->Eta())<0.7 && vec->Phi()>0 && vec->Phi()< phiMax) inEmc=kTRUE; | |
174 | vec = vec2; | |
175 | } | |
176 | hPartDiffPt->Fill(vec1->Pt()-vec2->Pt()); | |
177 | hPartDiffEta->Fill(vec1->Eta()-vec2->Eta()); | |
178 | dphi = vec1->Phi() - vec2->Phi(); | |
179 | dphi = fabs(TVector2::Phi_mpi_pi(dphi)); | |
180 | hPartDiffPhi->Fill(dphi); | |
181 | return inEmc; | |
182 | } | |
183 | ||
184 | void compareJetPartons(TVector3 *vecJet, TVector3 *vec1, TVector3 *vec2) | |
185 | {// 30-jan-2002 | |
186 | TVector3 *v, *vw; | |
187 | Double_t ang, angw; | |
188 | ||
189 | v = vec1; | |
190 | ang = vecJet->Angle((*v)); | |
191 | vw = vec2; | |
192 | angw = vecJet->Angle((*vw)); | |
193 | if(ang>angw) { | |
194 | ang = angw; | |
195 | v = vw; | |
196 | } | |
197 | hJetPartAng->Fill(ang); | |
198 | hDeta->Fill(vecJet->Eta() - v->Eta()); | |
199 | hDphi->Fill(vecJet->Phi() - v->Phi()); | |
200 | hDpt->Fill(vecJet->Pt() - v->Pt()); | |
201 | } | |
202 | ||
203 | void bookHist1() | |
204 | { | |
205 | gROOT->cd(); | |
206 | // Book histos | |
207 | hNjet= new TH1F("hNjet","# jets", 11, -0.5, 10.5); | |
208 | hJetEt = new TH1F("hJetEt","Energy of jet", 150, 0.0, 450.); | |
209 | hJetEta = new TH1F("hJetEta","#eta of jet", 180, -0.9, 0.9); | |
210 | hJetPhi = new TH1F("hJetPhi","#phi of jet", 62, -3.1, 3.1); | |
211 | hNtracksInJet = new TH1F("hNtracksInJet","n tracks", 30, 0.5, 29.5); | |
212 | hPtTracksInJet = new TH1F("hPtTracksInJet","p_{T} of tracks in jets cone", | |
213 | 100., 0., 100.); | |
214 | hTrackDR = new TH1F("hTrackDR","Track dR", 120., 0., 6.); | |
215 | #if !defined(__CINT__) | |
216 | lJets = _moveHistsToList("Hist for jets"); | |
217 | #endif | |
218 | // QA | |
219 | hJetPartAng = new TH1F("hJetPartAng","angle between jet and nearest parton", 100,0.0,0.2); | |
220 | hDeta = new TH1F("hDeta","#eta (jets) - #eta (partons)", 100, -.1, +.1); | |
221 | hDphi = new TH1F("hDphi","#phi (jets} - #phi (partons)", 100, -.1, +.1); | |
222 | hDpt = new TH1F("hDpt","E_{t} (jets) - E_{t} (partons)",100, -200., +200.); | |
223 | #if !defined(__CINT__) | |
224 | lCompPJ = _moveHistsToList("Hist for jet vs parton"); | |
225 | #endif | |
226 | } | |
227 | ||
228 | void bookHistPartonsAfterHardSc() | |
229 | { | |
230 | hPartPt = new TH2F("hPartPt","Partons p_{t}", 500, 30., 1030., 2, -0.5, 1.5); | |
231 | hPartEta = new TH2F("hPartEta","Partons #eta", 100, -5., 5., 2,-0.5,1.5); | |
232 | hPartPhi = new TH2F("hPartPhi","Partons #phi", 120,-0.2,2.2, 2,-0.5,1.5); | |
233 | hPartDiffPt = new TH1F("hPartDiffPt","p_{t} difference for partons",100, -100., 100.); | |
234 | hPartDiffEta = new TH1F("hPartDiffEta","#eta difference for partons",100, -5., 5.); | |
235 | hPartDiffPhi = new TH1F("hPartDiffPhi","#phi angle for partons",320, 0.0, 3.2); | |
236 | #if !defined(__CINT__) | |
237 | lPartons = _moveHistsToList("Hist for partons"); | |
238 | #endif | |
239 | } | |
240 | ||
241 | void drawPartons() | |
242 | { | |
243 | gStyle->SetOptStat(1111111); | |
244 | c3 = new TCanvas("c3","Partons canvas",400,10,600,700); | |
245 | c3->Divide(2,3); | |
246 | TH1D *hid; | |
247 | ||
248 | c3->cd(1); | |
249 | hid = hPartPt->ProjectionX("_pt1",1,1); _drawHist(hid, 2); | |
250 | hid = hPartPt->ProjectionX("_pt2",2,2); _drawHist(hid, 1, 2, "same"); | |
251 | c3->cd(2); _drawHist(hPartDiffPt,2); | |
252 | ||
253 | c3->cd(3); | |
254 | hid = hPartEta->ProjectionX("_eta1",1,1); _drawHist(hid, 2); | |
255 | hid = hPartEta->ProjectionX("_eta2",2,2); _drawHist(hid, 1, 2, "same"); | |
256 | c3->cd(4); _drawHist(hPartDiffEta,2); | |
257 | ||
258 | c3->cd(5); | |
259 | hid = hPartPhi->ProjectionX("_phi1",1,1); _drawHist(hid, 2); | |
260 | hid = hPartPhi->ProjectionX("_phi2",2,2); _drawHist(hid, 1, 2, "same"); | |
261 | c3->cd(6); _drawHist(hPartDiffPhi,2); | |
262 | } | |
263 | ||
264 | void fitGauss(TH1F* hid, Int_t opt) | |
265 | { // see /home/pavlinov/cosmic/pythia/geant/lin.C | |
266 | Double_t c = 2.5; | |
267 | Double_t xmax = hid->GetBinCenter(hid->GetMaximumBin()); | |
268 | Double_t rms = hid->GetRMS(); | |
269 | Double_t xmi=xmax -c*rms, xma=xmax+c*rms; | |
270 | _fit(111); | |
271 | if(!g1) { | |
272 | g1 = new TF1("g1","gaus", xmi, xma); | |
273 | } | |
274 | g1->SetLineColor(2); // red | |
275 | g1->SetLineWidth(3); | |
276 | ||
277 | TString stmp = fname(0,fname.Length()-5); // discard .root | |
278 | if(!c1) { | |
279 | c1 = new TCanvas("c1","Jet resolution",400,10,600,700); | |
280 | p1 = _pad(stmp.Data(),c1); | |
281 | } | |
282 | p1->cd(); p1->Clear(); | |
283 | hid->Fit("g1","R+"); | |
284 | ||
285 | if(opt) { | |
286 | _drawHist(hid, 2); | |
287 | if(opt>=10) { | |
288 | stmp += "_res.ps"; | |
289 | stmp.ReplaceAll("/FILE",""); | |
290 | p1->Print(stmp.Data()); | |
291 | } | |
292 | } | |
293 | } | |
294 | ||
295 | void cleanUpEvent() | |
296 | { | |
297 | /* ??? | |
298 | if(lPart[0]) { | |
299 | lPart[0]->Delete(); | |
300 | lPart[0]=0; | |
301 | lPart[1]->Delete(); | |
302 | lPart[1]=0; | |
303 | } | |
304 | */ | |
305 | } |