]>
Commit | Line | Data |
---|---|---|
a9e2aefa | 1 | TH1F *mass1 = new TH1F("mass1","Invariant Mass",120,0,12); |
2 | TH1F *mass2 = new TH1F("mass2","Invariant Mass",100,9.0,10.5); | |
3 | TH1F *mass3 = new TH1F("mass3","Invariant Mass",100,2.5,3.5); | |
4 | void MUONstraggling (Int_t evNumber1=0,Int_t evNumber2=0) | |
5 | { | |
6 | ///////////////////////////////////////////////////////////////////////// | |
7 | // This macro is a small example of a ROOT macro | |
8 | // illustrating how to read the output of GALICE | |
9 | // and do some analysis. | |
10 | // | |
11 | ///////////////////////////////////////////////////////////////////////// | |
12 | ||
13 | // Dynamically link some shared libs | |
14 | static Float_t xmuon, ymuon; | |
15 | ||
16 | if (gClassTable->GetID("AliRun") < 0) { | |
17 | gROOT->LoadMacro("loadlibs.C"); | |
18 | loadlibs(); | |
19 | } | |
20 | ||
21 | // Connect the Root Galice file containing Geometry, Kine and Hits | |
22 | ||
23 | TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); | |
24 | ||
25 | ||
26 | if (!file) { | |
27 | printf("\n Creating galice.root \n"); | |
28 | file = new TFile("galice.root"); | |
29 | } else { | |
30 | printf("\n galice.root found in file list"); | |
31 | } | |
32 | file->ls(); | |
33 | ||
34 | // Get AliRun object from file or create it if not on file | |
35 | if (!gAlice) { | |
36 | gAlice = (AliRun*)(file->Get("gAlice")); | |
37 | if (gAlice) printf("AliRun object found on file\n"); | |
38 | if (!gAlice) { | |
39 | printf("\n Create new gAlice object"); | |
40 | gAlice = new AliRun("gAlice","Alice test program"); | |
41 | } | |
42 | } | |
43 | ||
44 | ||
45 | // Create some histograms | |
46 | ||
47 | ||
48 | AliMUONChamber* iChamber; | |
49 | // | |
50 | // Loop over events | |
51 | // | |
52 | Int_t Nh=0; | |
53 | Int_t Nh1=0; | |
54 | for (Int_t nev=0; nev<= evNumber2; nev++) { | |
55 | Int_t nparticles = gAlice->GetEvent(nev); | |
56 | //cout << "nparticles " << nparticles <<endl; | |
57 | if (nev < evNumber1) continue; | |
58 | if (nparticles <= 0) return; | |
59 | ||
60 | AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); | |
61 | ||
62 | ||
63 | TTree *TH = gAlice->TreeH(); | |
64 | Int_t ntracks = TH->GetEntries(); | |
65 | // | |
66 | // Loop over tracks | |
67 | // | |
68 | ||
69 | Float_t pups[4]={0,0,0,0}; | |
70 | for (Int_t track=0; track<ntracks;track++) { | |
71 | Int_t Nelt=0; | |
72 | gAlice->ResetHits(); | |
73 | Int_t nbytes += TH->GetEvent(track); | |
74 | ||
75 | ||
76 | if (MUON) { | |
77 | for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); | |
78 | mHit; | |
79 | mHit=(AliMUONHit*)MUON->NextHit()) | |
80 | { | |
1e63a707 | 81 | Int_t nch = mHit->Chamber(); // chamber number |
82 | Float_t x = mHit->X(); // x-pos of hit | |
83 | Float_t y = mHit->Y(); // y-pos | |
84 | Float_t z = mHit->Z(); // y-pos | |
85 | Float_t p=mHit->Momentum(); | |
86 | Float_t px=mHit->Px(); | |
87 | Float_t py=mHit->Py(); | |
88 | Float_t pz=mHit->Pz(); | |
a9e2aefa | 89 | |
90 | if (nch != 1) continue; | |
91 | ||
1e63a707 | 92 | Int_t ipart = mHit->Particle(); |
a9e2aefa | 93 | TParticle *Part; |
1e63a707 | 94 | Int_t ftrack = mHit->Track(); |
95 | Part = gAlice->Particle(ftrack); | |
a9e2aefa | 96 | Int_t ipart = Part->GetPdgCode(); |
97 | TParticle *Mother; | |
98 | Float_t px0=Part->Px(); | |
99 | Float_t py0=Part->Py(); | |
100 | Float_t pz0=Part->Pz(); | |
101 | Float_t e0=Part->Energy(); | |
102 | ||
103 | if (ipart == kMuonPlus || ipart == kMuonMinus) { | |
104 | // | |
105 | // Branson Correction | |
106 | // | |
107 | Float_t zch=505.; | |
108 | Float_t r=TMath::Sqrt(x*x+y*y); | |
109 | Float_t zb; | |
110 | ||
111 | if (r<26.3611) { | |
112 | zb=466.; | |
113 | } else { | |
114 | zb=441.; | |
115 | } | |
116 | ||
117 | Float_t xb=x-(zch-zb)*px/pz; | |
118 | Float_t yb=y-(zch-zb)*py/pz; | |
119 | Float_t tx=xb/zb; | |
120 | Float_t ty=yb/zb; | |
121 | pz=zb*p/TMath::Sqrt(zb*zb+xb*xb+yb*yb); | |
122 | px=pz*tx; | |
123 | py=pz*ty; | |
124 | ||
125 | // | |
126 | // Energy Correction | |
127 | // | |
128 | // | |
129 | Float_t corr=(p+CorrectP(p,mHit->fTheta))/p; | |
130 | pups[0]+=p*corr; | |
131 | pups[1]+=px*corr; | |
132 | pups[2]+=py*corr; | |
133 | pups[3]+=pz*corr; | |
134 | } | |
135 | } // hits | |
136 | } // if MUON | |
137 | } // tracks | |
138 | Float_t mass=TMath::Sqrt(pups[0]*pups[0]-pups[1]*pups[1]-pups[2]*pups[2]-pups[3]*pups[3]); | |
139 | mass1->Fill(mass, 1.); | |
140 | mass2->Fill(mass, 1.); | |
141 | mass3->Fill(mass, 1.); | |
142 | } // event | |
143 | //Create a canvas, set the view range, show histograms | |
144 | TCanvas *c1 = new TCanvas("c1","Vertices from electrons and positrons",400,10,600,700); | |
145 | c1->Divide(2,2); | |
146 | ||
147 | c1->cd(1); | |
148 | mass1->SetFillColor(42); | |
149 | mass1->SetXTitle("Mass (GeV)"); | |
150 | mass1->Draw(); | |
151 | ||
152 | c1->cd(2); | |
153 | mass2->SetFillColor(42); | |
154 | mass2->SetXTitle("Mass (GeV)"); | |
155 | mass2->Draw(); | |
156 | ||
157 | c1->cd(3); | |
158 | mass3->SetFillColor(42); | |
159 | mass3->SetXTitle("Mass (GeV)"); | |
160 | mass3->Draw(); | |
161 | ||
162 | } | |
163 | ||
164 | ||
165 | ||
166 | Float_t CorrectP(Float_t p, Float_t theta) | |
167 | { | |
168 | if (theta<3.) { | |
169 | //W | |
170 | if (p<15) { | |
171 | return 2.737+0.0494*p-0.001123*p*p; | |
172 | } else { | |
173 | return 3.0643+0.01346*p; | |
174 | } | |
175 | } else { | |
176 | //Pb | |
177 | if (p<15) { | |
178 | return 2.1380+0.0351*p-0.000853*p*p; | |
179 | } else { | |
180 | return 2.407+0.00702*p; | |
181 | } | |
182 | } | |
183 | } | |
184 | ||
185 | ||
186 | ||
187 | ||
188 | ||
189 |