Energy correction for default absorber configuration.
[u/mrichter/AliRoot.git] / MUON / MUONstraggling.C
CommitLineData
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);
4void 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
158Float_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