]>
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 | } | |
a9e2aefa | 32 | // Get AliRun object from file or create it if not on file |
33 | if (!gAlice) { | |
34 | gAlice = (AliRun*)(file->Get("gAlice")); | |
35 | if (gAlice) printf("AliRun object found on file\n"); | |
36 | if (!gAlice) { | |
37 | printf("\n Create new gAlice object"); | |
38 | gAlice = new AliRun("gAlice","Alice test program"); | |
39 | } | |
40 | } | |
41 | ||
42 | ||
43 | // Create some histograms | |
44 | ||
45 | ||
46 | AliMUONChamber* iChamber; | |
47 | // | |
48 | // Loop over events | |
49 | // | |
50 | Int_t Nh=0; | |
51 | Int_t Nh1=0; | |
52 | for (Int_t nev=0; nev<= evNumber2; nev++) { | |
53 | Int_t nparticles = gAlice->GetEvent(nev); | |
54 | //cout << "nparticles " << nparticles <<endl; | |
55 | if (nev < evNumber1) continue; | |
56 | if (nparticles <= 0) return; | |
57 | ||
58 | AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON"); | |
59 | ||
60 | ||
61 | TTree *TH = gAlice->TreeH(); | |
62 | Int_t ntracks = TH->GetEntries(); | |
63 | // | |
64 | // Loop over tracks | |
65 | // | |
66 | ||
67 | Float_t pups[4]={0,0,0,0}; | |
68 | for (Int_t track=0; track<ntracks;track++) { | |
69 | Int_t Nelt=0; | |
70 | gAlice->ResetHits(); | |
71 | Int_t nbytes += TH->GetEvent(track); | |
72 | ||
73 | ||
74 | if (MUON) { | |
75 | for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); | |
76 | mHit; | |
77 | mHit=(AliMUONHit*)MUON->NextHit()) | |
78 | { | |
1e63a707 | 79 | Int_t nch = mHit->Chamber(); // chamber number |
80 | Float_t x = mHit->X(); // x-pos of hit | |
81 | Float_t y = mHit->Y(); // y-pos | |
82 | Float_t z = mHit->Z(); // y-pos | |
9a773b96 | 83 | Float_t p = mHit->Momentum(); |
84 | Float_t px = mHit->Px(); | |
85 | Float_t py = mHit->Py(); | |
86 | Float_t pz = mHit->Pz(); | |
a9e2aefa | 87 | |
88 | if (nch != 1) continue; | |
89 | ||
1e63a707 | 90 | Int_t ipart = mHit->Particle(); |
a9e2aefa | 91 | TParticle *Part; |
1e63a707 | 92 | Int_t ftrack = mHit->Track(); |
93 | Part = gAlice->Particle(ftrack); | |
a9e2aefa | 94 | Int_t ipart = Part->GetPdgCode(); |
95 | TParticle *Mother; | |
9a773b96 | 96 | Float_t px0 = Part->Px(); |
97 | Float_t py0 = Part->Py(); | |
98 | Float_t pz0 = Part->Pz(); | |
99 | Float_t p0 = Part->P(); | |
100 | ||
a9e2aefa | 101 | Float_t e0=Part->Energy(); |
102 | ||
103 | if (ipart == kMuonPlus || ipart == kMuonMinus) { | |
104 | // | |
105 | // Branson Correction | |
106 | // | |
9a773b96 | 107 | Float_t zch = z; |
108 | Float_t r = TMath::Sqrt(x*x+y*y); | |
a9e2aefa | 109 | Float_t zb; |
110 | ||
111 | if (r<26.3611) { | |
9a773b96 | 112 | zb = 466.; |
a9e2aefa | 113 | } else { |
9a773b96 | 114 | zb = 441.; |
a9e2aefa | 115 | } |
9a773b96 | 116 | |
117 | Float_t xb = x-(zch-zb)*px/pz; | |
118 | Float_t yb = y-(zch-zb)*py/pz; | |
119 | Float_t pzB = p * zb/TMath::Sqrt(zb*zb+yb*yb+xb*xb); | |
120 | Float_t pxB = pzB * xb/zb; | |
121 | Float_t pyB = pzB * yb/zb; | |
122 | Float_t theta | |
123 | = TMath::ATan2(TMath::Sqrt(pxB*pxB+pyB*pyB), pzB) * | |
124 | 180./TMath::Pi(); | |
a9e2aefa | 125 | |
a9e2aefa | 126 | |
127 | // | |
128 | // Energy Correction | |
129 | // | |
130 | // | |
9a773b96 | 131 | |
132 | Float_t corr=(p+CorrectP(p, theta))/p; | |
a9e2aefa | 133 | pups[0]+=p*corr; |
9a773b96 | 134 | pups[1]+=pxB*corr; |
135 | pups[2]+=pyB*corr; | |
136 | pups[3]+=pzB*corr; | |
a9e2aefa | 137 | } |
138 | } // hits | |
139 | } // if MUON | |
140 | } // tracks | |
9a773b96 | 141 | Float_t mass = |
142 | TMath::Sqrt(pups[0]*pups[0] | |
143 | -pups[1]*pups[1]-pups[2]*pups[2]-pups[3]*pups[3]); | |
a9e2aefa | 144 | mass1->Fill(mass, 1.); |
145 | mass2->Fill(mass, 1.); | |
146 | mass3->Fill(mass, 1.); | |
147 | } // event | |
148 | //Create a canvas, set the view range, show histograms | |
9a773b96 | 149 | TCanvas *c1 = new TCanvas("c1","Mass Spectra",400,10,600,700); |
150 | ||
a9e2aefa | 151 | mass2->SetFillColor(42); |
152 | mass2->SetXTitle("Mass (GeV)"); | |
153 | mass2->Draw(); | |
a9e2aefa | 154 | } |
155 | ||
156 | ||
157 | ||
158 | Float_t CorrectP(Float_t p, Float_t theta) | |
159 | { | |
9a773b96 | 160 | Float_t c; |
161 | ||
a9e2aefa | 162 | if (theta<3.) { |
163 | //W | |
164 | if (p<15) { | |
9a773b96 | 165 | c = 2.916+0.0156*p-0.00000866*p*p; |
a9e2aefa | 166 | } else { |
9a773b96 | 167 | c = 2.998+0.0137*p; |
a9e2aefa | 168 | } |
169 | } else { | |
9a773b96 | 170 | //Cu+Fe+Cc+C |
a9e2aefa | 171 | if (p<15) { |
9a773b96 | 172 | c = 2.464+0.00877*p-0.0000106*p*p; |
a9e2aefa | 173 | } else { |
9a773b96 | 174 | c = 2.496+0.00817*p; |
a9e2aefa | 175 | } |
176 | } | |
9a773b96 | 177 | // |
178 | // | |
179 | c*=0.95; | |
180 | // | |
181 | return c; | |
a9e2aefa | 182 | } |
183 | ||
184 | ||
185 | ||
186 | ||
187 | ||
188 |