- Welding section on absorber side (LHCVC2C_001)
[u/mrichter/AliRoot.git] / MUON / MUONacc.C
CommitLineData
9efe4ffa 1Bool_t SelectThis(AliMUONHit* mHit);
2
3Bool_t SelectThis(Float_t x, Float_t y, Float_t d)
4{
04611e14 5// Float_t x = mHit->X(); // x-pos
6// Float_t y = mHit->Y(); // y-pos
9efe4ffa 7 if (TMath::Abs(x) > d/2. && TMath::Abs(y) > d/2.) {
8 return 1;
9 } else {
10 return 0;
11 }
12}
13
04611e14 14Bool_t SelectThat(Float_t x, Float_t y, Float_t d)
15{
16 const Float_t phi0=60.*TMath::Pi()/180.;
17
18 if (y<0.) return 1;
19 if (x<0.) x=-x;
20
21 Float_t phi = TMath::ATan2(y,x);
22 Float_t dist = TMath::Sin(phi0-phi)*TMath::Sqrt(x*x+y*y);
23 if (TMath::Abs(dist) > d) {
24 return 1;
25 } else {
26 return 0 ;
27 }
28}
9efe4ffa 29
30void MUONacc (Float_t d=3., Int_t evNumber1=0, Int_t evNumber2=0)
31{
32// Dynamically link some shared libs
33 Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
34 56.1, 56.1, 59.6, 59.6};
35 Float_t rmax[14] = {91.5, 91.5, 122.5, 122.5, 158.3, 158.3, 260.0,
36 260.0, 260.0, 260.0, 293., 293, 293., 293.};
37
38
39 if (gClassTable->GetID("AliRun") < 0) {
40 gROOT->LoadMacro("loadlibs.C");
41 loadlibs();
42 }
43// Connect the Root Galice file containing Geometry, Kine and Hits
44
45 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
46
47
48 if (!file) {
49 printf("\n Creating galice.root \n");
50 file = new TFile("galice.root");
51 } else {
52 printf("\n galice.root found in file list");
53 }
54 file->ls();
55//
56 TDatabasePDG* DataBase = new TDatabasePDG();
57
58// Get AliRun object from file or create it if not on file
59 if (!gAlice) {
60 gAlice = (AliRun*)(file->Get("gAlice"));
61 if (gAlice) printf("AliRun object found on file\n");
62 if (!gAlice) {
63 printf("\n create new gAlice object");
64 gAlice = new AliRun("gAlice","Alice test program");
65 }
66 }
67
68// Create some histograms
69
70 TH1F *theta = new TH1F("theta","Theta distribution",180,0,180);
71 TH1F *emult = new TH1F("emult","Event Multiplicity",100,0,1000);
04611e14 72 TH2F *supp = new TH2F("supp", "Diselected",100, -200., 200, 100, -200., 200.);
9efe4ffa 73
74 AliMUONChamber* iChamber;
75 AliSegmentation* seg;
76
77//
78// Loop over events
79//
80
81 for (Int_t nev=0; nev<= evNumber2; nev++) {
82 Int_t nparticles = gAlice->GetEvent(nev);
83 if (nev < evNumber1) continue;
84 if (nparticles <= 0) return;
85
86 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
87
88
89 TTree *TH = gAlice->TreeH();
90 Int_t ntracks = TH->GetEntries();
91 TClonesArray *fPartArray = gAlice->Particles();
92 Int_t npart = fPartArray->GetEntriesFast();
93
94 Int_t EvMult=0;
95 Int_t Nsel=0;
96//
97// Loop over primary particles (jpsi. upsilon, ...)
98//
99 Int_t prim=0;
100
04611e14 101 for (Int_t part=0; part<30000; part+=3) {
102 TParticle *MPart = gAlice->Particle(part);
103 Int_t mpart = MPart->GetPdgCode();
9efe4ffa 104 Int_t child1 = MPart->GetFirstDaughter();
105 Int_t child2 = MPart->GetLastDaughter();
106 Int_t mother = MPart->GetFirstMother();
107//
108// Loop over muons
109//
110 Int_t muon1=2*prim;
111 Int_t muon2=muon1+1;
112
113// printf("\n PDG Code %d %d %d %d %d %d %d %d \n",
114// part, mpart, child1, child2, mother,prim, muon1,muon2);
115
116
117 Bool_t selected[2]={1,1};
118 Int_t muons=0;
119 for (Int_t track=muon1; track<=muon2;track++) {
120 gAlice->ResetHits();
121 Int_t nbytes += TH->GetEvent(track);
122
123 if (MUON) {
124//
125// Loop over hits
126//
127 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
128 mHit;
129 mHit=(AliMUONHit*)MUON->NextHit())
130 {
04611e14 131 Int_t nch = mHit->Chamber(); // chamber number
132 Float_t x = mHit->X(); // x-pos of hit
133 Float_t y = mHit->Y(); // y-pos
134 Float_t z = mHit->Z(); // z-pos
135 Float_t Eloss = mHit->Eloss(); // energy loss
136 Float_t Theta = mHit->Theta(); // theta
137 Float_t Particle = mHit->Particle(); // Particle type
138 Int_t itrack = Int_t(mHit->Track());
139
140 TParticle *thePart = gAlice->Particle(itrack);
9efe4ffa 141 Float_t pTheta=thePart->Theta();
142
143 if (Particle != kMuonPlus && Particle != kMuonMinus) continue;
144
1e63a707 145 Float_t P = mhit->Momentum();
9efe4ffa 146 Float_t R = TMath::Sqrt(x*x+y*y);
147 TParticlePDG* Part = DataBase->GetParticle(Particle);
148 Double_t mass = Part->Mass();
149 theta->Fill(pTheta*180./TMath::Pi(),1.);
150
151 if (nch>14) continue;
152// printf("\n %f %f %f %d\n ", R, rmin[nch-1], rmax[nch-1], nch);
153// Geometrical acceptance of frames
154// selected[muons] = selected[muons]&&SelectThis(x,y,d);
155// Cut on default rmin and rmax
156// selected[muons] = selected[muons] &&( R > rmin[nch-1] && R < rmax[nch-1]);
157// Cut on rmin
158// selected[muons] = selected[muons] &&( R > rmin[nch-1]);
159// Cut on rmax
160// selected[muons] = selected[muons] &&( R < rmax[nch-1]);
161// selected[muons] = selected[muons] &&( R > z*TMath::Tan(2.0*TMath::Pi()/180.));
162// Bool_t ok = (z<970 && ( R > z*TMath::Tan(2.0*TMath::Pi()/180.))) ||
163// (z>970 && ( R > 970 * TMath::Tan(2.0*TMath::Pi()/180.)
164// +(z-970)* TMath::Tan(1.6*TMath::Pi()/180.)));
165//
04611e14 166// Bool_t ok = (z<970 && ( R < z*TMath::Tan(9.0*TMath::Pi()/180.))) ||
167// (z>970 && ( R < 970 * TMath::Tan(9.0*TMath::Pi()/180.)
168// +(z-970)* TMath::Tan(12.0*TMath::Pi()/180.)));
169 if (nch == 4 || nch == 7) {
170 Bool_t ok = SelectThat(x,y,d);
171 } else {
172 Bool_t ok = kTRUE;
173 }
174
175
176 selected[muons] = selected[muons] && ok;
177 if (!ok) {
178 supp->Fill(x,y,1.);
179
180 }
9efe4ffa 181
9efe4ffa 182
183 } // hit loop
184 } // if MUON
185 muons++;
186 } // track loop
187 if (selected[0] && selected[1]) Nsel++;
188 prim++;
189 } // primary loop
190 printf("\n Selected %d %d \n", Nsel, prim);
191 }
192
193//Create a canvas, set the view range, show histograms
194 Int_t k;
195 TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
04611e14 196 c1->Divide(2,2);
9efe4ffa 197 c1->cd(1);
198 theta->Draw();
04611e14 199 c1->cd(2);
200 supp->Draw();
9efe4ffa 201
202}
203
204
205
206
207
208
209
210
211