Code from MUON-dev joined
[u/mrichter/AliRoot.git] / MUON / MUONtestabso.C
CommitLineData
a897a37a 1void 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;
140 Float_t thetax=TMath::ATan2(mHit->fCxHit, mHit->fPTot);
141 Float_t thetay=TMath::ATan2(mHit->fCyHit, mHit->fPTot);
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
253Float_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