]>
Commit | Line | Data |
---|---|---|
a897a37a | 1 | void MUONtestabso (Int_t evNumber1=0,Int_t evNumber2=0) |
2 | { | |
3 | ///////////////////////////////////////////////////////////////////////// | |
4 | // This macro is a small example of a ROOT macro | |
5 | // illustrating how to read the output of GALICE | |
6 | // and do some analysis. | |
7 | // | |
8 | ///////////////////////////////////////////////////////////////////////// | |
9 | ||
10 | // Dynamically link some shared libs | |
a9e2aefa | 11 | static Float_t xmuon, ymuon; |
12 | ||
a897a37a | 13 | if (gClassTable->GetID("AliRun") < 0) { |
14 | gROOT->LoadMacro("loadlibs.C"); | |
15 | loadlibs(); | |
16 | } | |
17 | ||
18 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
19 | ||
20 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
21 | ||
22 | ||
23 | if (!file) { | |
24 | printf("\n Creating galice.root \n"); | |
25 | file = new TFile("galice.root"); | |
26 | } else { | |
27 | printf("\n galice.root found in file list"); | |
28 | } | |
29 | file->ls(); | |
30 | ||
31 | // Get AliRun object from file or create it if not on file | |
32 | if (!gAlice) { | |
33 | gAlice = (AliRun*)(file->Get("gAlice")); | |
34 | if (gAlice) printf("AliRun object found on file\n"); | |
35 | if (!gAlice) { | |
a9e2aefa | 36 | printf("\n Create new gAlice object"); |
a897a37a | 37 | gAlice = new AliRun("gAlice","Alice test program"); |
38 | } | |
39 | } | |
40 | ||
a897a37a | 41 | |
42 | // Create some histograms | |
43 | ||
44 | ||
a9e2aefa | 45 | TH1F *zv = new TH1F("zv","z vertex" ,140,470.,505.); |
46 | TH1F *xv = new TH1F("xv","x vertex" ,140,-70.,70.); | |
47 | TH1F *yv = new TH1F("yv","y vertex", 140,-70.,70.); | |
48 | ||
49 | TH1F *ip = new TH1F("ip","geant part",50,0.,10.); | |
50 | TH2F *rzv = new TH2F("rzv","R-z vert",100,502,504.,100,0.,100.); | |
51 | ||
52 | TH1F *ptUps = new TH1F("ptUps","pT Upsilon",50,0.,25.); | |
53 | TH1F *hde = new TH1F("hde","dE",200,0.,50); | |
54 | TH1F *hde2 = new TH1F("hde2","dE",100,1.5,5.5); | |
55 | TH2F *hdevsn = new TH2F("hdevsn","dE vs N electron",100,0.,15., 20, 0.,20.); | |
a897a37a | 56 | |
a9e2aefa | 57 | TH1F *ekine = new TH1F("ekine","E_kin electrons",70,-5,2); |
58 | TH1F *etheta = new TH1F("etheta","Theta electrons",90,0,90); | |
59 | TH1F *edr = new TH1F("edr","Distance to muon",100,0,10); | |
a897a37a | 60 | |
a9e2aefa | 61 | TH1F *de = new TH1F("de","correction",100,-1,1); |
62 | ||
63 | TH1F *dtheta = new TH1F("dtheta","Delta Theta" ,200,-5.,5.); | |
64 | AliMUONChamber* iChamber; | |
a897a37a | 65 | // |
66 | // Loop over events | |
67 | // | |
68 | Int_t Nh=0; | |
69 | Int_t Nh1=0; | |
a9e2aefa | 70 | Int_t Nel1=0; |
71 | Int_t Nel2=0; | |
72 | Int_t Nel3=0; | |
73 | Int_t Nel4=0; | |
74 | ||
a897a37a | 75 | for (Int_t nev=0; nev<= evNumber2; nev++) { |
a9e2aefa | 76 | //cout << "nev " << nev <<endl; |
a897a37a | 77 | Int_t nparticles = gAlice->GetEvent(nev); |
a9e2aefa | 78 | //cout << "nparticles " << nparticles <<endl; |
a897a37a | 79 | if (nev < evNumber1) continue; |
80 | if (nparticles <= 0) return; | |
a9e2aefa | 81 | |
a897a37a | 82 | AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); |
a9e2aefa | 83 | |
a897a37a | 84 | |
85 | TTree *TH = gAlice->TreeH(); | |
86 | Int_t ntracks = TH->GetEntries(); | |
87 | // | |
88 | // Loop over tracks | |
89 | // | |
a9e2aefa | 90 | |
91 | Float_t dE; | |
92 | ||
a897a37a | 93 | for (Int_t track=0; track<ntracks;track++) { |
a9e2aefa | 94 | Int_t Nelt=0; |
95 | printf("\n nev %d %d %d %d %d\n ", nev, Nel1, Nel2, Nel3, Nel4); | |
a897a37a | 96 | gAlice->ResetHits(); |
97 | Int_t nbytes += TH->GetEvent(track); | |
98 | if (MUON) { | |
a9e2aefa | 99 | for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); |
a897a37a | 100 | mHit; |
a9e2aefa | 101 | mHit=(AliMUONHit*)MUON->NextHit()) |
a897a37a | 102 | { |
103 | Int_t nch = mHit->fChamber; // chamber number | |
104 | Float_t x = mHit->fX; // x-pos of hit | |
105 | Float_t y = mHit->fY; // y-pos | |
106 | Float_t z = mHit->fZ; // y-pos | |
a9e2aefa | 107 | |
108 | if (nch != 1) continue; | |
109 | ||
110 | Int_t ipart = mHit->fParticle; | |
111 | TClonesArray *fPartArray = gAlice->Particles(); | |
112 | TParticle *Part; | |
113 | Int_t ftrack = mHit->fTrack; | |
114 | Part = (TParticle*) fPartArray->UncheckedAt(ftrack); | |
115 | Int_t ipart = Part->GetPdgCode(); | |
116 | ip->Fill((float)ipart); | |
117 | TParticle *Mother; | |
118 | ||
119 | Float_t px0=Part->Px(); | |
120 | Float_t py0=Part->Py(); | |
121 | Float_t pz0=Part->Pz(); | |
122 | Float_t thetax0=TMath::ATan2(px0,pz0); | |
123 | Float_t thetay0=TMath::ATan2(py0,pz0); | |
124 | ||
125 | if (ipart == kMuonPlus || ipart == kMuonMinus) { | |
126 | // Int_t imo = Part->GetFirstMother(); | |
127 | // Mother = (TParticle*) fPartArray->UncheckedAt(imo); | |
128 | // Float_t pt = Mother->Pt(); | |
129 | // ptUps->Fill(pt, (float) 1); | |
130 | Float_t E=Part->Energy(); | |
131 | Float_t Eloc=mHit->fPTot; | |
132 | Float_t corr=Eloc+CorrectP(Eloc,mHit->fTheta); | |
133 | printf("\n %f %f %f", E, Eloc, corr); | |
134 | de->Fill(E-corr,1.); | |
135 | dE = E-Eloc; | |
136 | if (dE<50) hde->Fill(dE, (float) 1); | |
137 | if (dE<5.5) hde2->Fill(dE, (float) 1); | |
138 | xmuon=mHit->fX; | |
139 | ymuon=mHit->fY; | |
1e63a707 | 140 | Float_t thetax=TMath::ATan2(mHit->Px(), mHit->Momentum()); |
141 | Float_t thetay=TMath::ATan2(mHit->Py(), mHit->Momentum()); | |
a9e2aefa | 142 | dtheta->Fill((thetax-thetax0)*1000., 1.); |
143 | dtheta->Fill((thetay-thetay0)*1000., 1.); | |
144 | } | |
145 | ||
146 | if (ipart == kElectron || ipart == kPositron) { | |
147 | ||
148 | ||
149 | Float_t xvert = Part->Vx(); // vertex | |
150 | Float_t yvert = Part->Vy(); // vertex | |
151 | Float_t zvert = Part->Vz(); // z vertex | |
152 | if (zvert < 503 && mHit->fTheta<90) { | |
153 | Nelt++; | |
154 | Float_t px = Part->Px(); | |
155 | Float_t py = Part->Py(); | |
156 | Float_t pz = Part->Pz(); | |
157 | Float_t Ek = Part->Energy()-Part->GetMass(); | |
158 | ||
159 | Int_t imo = Part->GetFirstMother(); | |
160 | Mother = (TParticle*) fPartArray->UncheckedAt(imo); | |
161 | Int_t imot = Mother->GetPdgCode(); | |
162 | xv->Fill(xvert); | |
163 | yv->Fill(yvert); | |
164 | zv->Fill(zvert); | |
165 | Float_t rvert=TMath::Sqrt(xvert*xvert+yvert*yvert); | |
166 | rzv->Fill(zvert,rvert); | |
167 | ekine->Fill(TMath::Log10(Ek),1.); | |
168 | etheta->Fill(mHit->fTheta,1.); | |
169 | Float_t ex=mHit->fX; | |
170 | Float_t ey=mHit->fY; | |
171 | dr=TMath::Sqrt((ex-xmuon)*(ex-xmuon)+(ey-ymuon)*(ey-ymuon)); | |
172 | edr->Fill(dr,1.); | |
173 | ||
174 | } | |
175 | } | |
176 | ||
177 | } // hits | |
178 | } // if MUON | |
179 | if (Nelt == 1) Nel1++; | |
180 | if (Nelt == 2) Nel2++; | |
181 | if (Nelt == 3) Nel3++; | |
182 | if (Nelt > 3) Nel4++; | |
183 | hdevsn->Fill(dE, (float) Nelt, (float) 1); | |
184 | ||
185 | } // tracks | |
186 | ||
187 | } // event | |
a897a37a | 188 | |
189 | ||
190 | //Create a canvas, set the view range, show histograms | |
191 | TCanvas *c1 = new TCanvas("c1","Vetices from electrons and positrons",400,10,600,700); | |
a9e2aefa | 192 | c1->Divide(2,2); |
a897a37a | 193 | |
a9e2aefa | 194 | c1->cd(1); |
a897a37a | 195 | ip->SetFillColor(42); |
196 | ip->SetXTitle("ipart"); | |
197 | ip->Draw(); | |
198 | ||
a9e2aefa | 199 | c1->cd(2); |
a897a37a | 200 | xv->SetFillColor(42); |
201 | xv->SetXTitle("xvert"); | |
202 | xv->Draw(); | |
203 | ||
a9e2aefa | 204 | c1->cd(3); |
a897a37a | 205 | yv->SetFillColor(42); |
206 | yv->SetXTitle("yvert"); | |
207 | yv->Draw(); | |
208 | ||
a9e2aefa | 209 | c1->cd(4); |
a897a37a | 210 | zv->SetFillColor(42); |
211 | zv->SetXTitle("zvert"); | |
212 | zv->Draw(); | |
213 | ||
a9e2aefa | 214 | TCanvas *c2 = new TCanvas("c2"," ",400,10,600,700); |
215 | c2->Divide(2,2); | |
216 | c2->cd(1); | |
a897a37a | 217 | rzv->SetXTitle("zvert"); |
218 | rzv->SetYTitle("rvert"); | |
219 | rzv->Draw(); | |
a9e2aefa | 220 | |
221 | c2->cd(2); | |
222 | ptUps->SetXTitle("pt"); | |
223 | ptUps->Draw(); | |
224 | ||
225 | c2->cd(3); | |
226 | hde->SetXTitle("dE"); | |
227 | hde->Draw(); | |
228 | ||
229 | ||
230 | c2->cd(4); | |
231 | hde2->SetXTitle("dE"); | |
232 | hde2->Draw(); | |
233 | TCanvas *c3 = new TCanvas("c3"," ",400,10,600,700); | |
234 | c3->Divide(2,2); | |
235 | c3->cd(1); | |
236 | ekine->SetXTitle("E_kin"); | |
237 | ekine->Draw(); | |
238 | ||
239 | c3->cd(2); | |
240 | etheta->SetXTitle("Theta"); | |
241 | etheta->Draw(); | |
242 | ||
243 | c3->cd(3); | |
244 | edr->SetXTitle("Distance to muon"); | |
245 | edr->Draw(); | |
246 | ||
247 | c3->cd(4); | |
248 | dtheta->SetXTitle(" "); | |
249 | dtheta->Draw(); | |
250 | ||
251 | } | |
252 | ||
253 | Float_t CorrectP(Float_t p, Float_t theta) | |
254 | { | |
255 | printf("\n %f%", theta); | |
256 | ||
257 | if (theta<3.) { | |
258 | //W | |
259 | if (p<15) { | |
260 | return 2.737+0.0494*p-0.001123*p*p; | |
261 | } else { | |
262 | return 3.0643+0.01346*p; | |
263 | } | |
264 | } else { | |
265 | //Pb | |
266 | if (p<15) { | |
267 | return 2.1380+0.0351*p-0.000853*p*p; | |
268 | } else { | |
269 | return 2.407+0.00702*p; | |
270 | } | |
271 | } | |
a897a37a | 272 | } |
a9e2aefa | 273 |