Code from MUON-dev joined
[u/mrichter/AliRoot.git] / MUON / MUONhits.C
1 void MUONhits (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ic=1) 
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                    
11     Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
12                         56.1, 56.1, 59.6, 59.6};
13     Float_t rmax[14] = {81.6, 81.6, 109.3, 109.3, 154.4, 154.4, 197.8, 
14                         197.8, 229.5, 229.5, 254., 254, 270., 270.};
15     
16
17     if (gClassTable->GetID("AliRun") < 0) {
18         gROOT->LoadMacro("loadlibs.C");
19         loadlibs();
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     TDatabasePDG*  DataBase = new TDatabasePDG();
35     
36 // Get AliRun object from file or create it if not on file
37     if (!gAlice) {
38         gAlice = (AliRun*)(file->Get("gAlice"));
39         if (gAlice) printf("AliRun object found on file\n");
40         if (!gAlice) {
41             printf("\n create new gAlice object");
42             gAlice = new AliRun("gAlice","Alice test program");
43         }
44     }
45     
46 //  Create some histograms
47
48     TH1F *Hits[14],   *HitDensity[14];
49     TH1F *Hits_g[14], *HitDensity_g[14];
50     TH1F *Hits_n[14], *HitDensity_n[14];
51     TH1F *Mom[14], *Mom_g[14], *Mom_n[14];
52     
53     for (Int_t i=0; i<14; i++) {
54         Hits[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
55         Hits[i]->SetFillColor(0);
56         Hits[i]->SetXTitle("R (cm)");
57         
58         HitDensity[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
59         HitDensity[i]->SetFillColor(0);
60         HitDensity[i]->SetXTitle("R (cm)");
61
62         Hits_g[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
63         Hits_g[i]->SetFillColor(0);
64         Hits_g[i]->SetXTitle("R (cm)");
65
66         HitDensity_g[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
67         HitDensity_g[i]->SetFillColor(0);
68         HitDensity_g[i]->SetXTitle("R (cm)");
69
70         Mom[i]  =  new TH1F("mom","Energy",70,-6,1);
71         Mom[i]->SetFillColor(0);
72         Mom[i]->SetXTitle("E (GeV)");
73
74         Mom_g[i]  =  new TH1F("mom","Energy",70,-6,1);
75         Mom_g[i]->SetFillColor(0);
76         Mom_g[i]->SetXTitle("E (GeV)");
77
78         Mom_n[i]  =  new TH1F("mom","Energy",70,-6,1);
79         Mom_n[i]->SetFillColor(0);
80         Mom_n[i]->SetXTitle("E (GeV)");
81
82         Hits_n[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
83         Hits_n[i]->SetFillColor(0);
84         Hits_n[i]->SetXTitle("R (cm)");
85
86         HitDensity_n[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
87         HitDensity_n[i]->SetFillColor(0);
88         HitDensity_n[i]->SetXTitle("R (cm)");
89     }
90
91     TH1F *theta   =  new TH1F("theta","Theta distribution",180,0,180);
92  
93     TH1F *emult   =  new TH1F("emult","Event Multiplicity",100,0,1000);   
94     AliMUONChamber*  iChamber;
95     AliMUONSegmentation* seg;
96     
97 //
98 //   Loop over events 
99 //
100     Float_t dpx[2], dpy[2];
101     Int_t mdig =0;
102     Int_t nhit=0;
103     
104     for (Int_t nev=0; nev<= evNumber2; nev++) {
105         Int_t nparticles = gAlice->GetEvent(nev);
106         if (nev < evNumber1) continue;
107         if (nparticles <= 0) return;
108         
109         AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
110         printf("\n track %d %d \n ", nev, MUON);
111         
112         TTree *TH = gAlice->TreeH();
113         Int_t ntracks = TH->GetEntries();
114 //
115 //   Loop over events
116 //
117         Int_t EvMult=0;
118        
119         for (Int_t track=0; track<ntracks;track++) {
120             gAlice->ResetHits();
121             Int_t nbytes += TH->GetEvent(track);
122             if (MUON)  {
123 //
124 //  Loop over hits
125 //
126                 Int_t NperCh=0;
127                
128                 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
129                     mHit;
130                     mHit=(AliMUONHit*)MUON->NextHit()) 
131                 {
132                     Int_t   nch   =    mHit->fChamber;  // chamber number
133                     Float_t x     =    mHit->fX;        // x-pos of hit
134                     Float_t y     =    mHit->fY;        // y-pos
135                     Float_t Eloss =    mHit->fEloss;
136                     Float_t Theta =    mHit->fTheta;
137                     Float_t Particle = mHit->fParticle;
138                     Float_t P    =     
139                         TMath::Sqrt(mHit->fCxHit*mHit->fCxHit+
140                                     mHit->fCyHit*mHit->fCyHit+
141                                     mHit->fCzHit*mHit->fCzHit);
142                     TParticlePDG* Part = DataBase->GetParticle(Particle);
143                     Double_t mass = Part->Mass();
144                     
145                     if (nch >13) continue;
146                     if (nch ==1) EvMult++;
147                     if (mHit->fAge > 5.e-6) continue;
148                    
149                     Float_t r=TMath::Sqrt(x*x+y*y);
150                     if (nch ==0) continue;
151                     
152                     if (r < rmin[nch-1]) continue;
153                     if (r > rmax[nch-1]) continue;
154                    
155                     Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);
156
157                     Float_t Ekin=TMath::Sqrt(P*P+mass*mass)-mass;
158                     if (Particle == 2112) {
159                         Hits_n[nch-1]->Fill(r,(float) 1);
160                         HitDensity_n[nch-1]->Fill(r,wgt);
161                         Mom_n[nch-1]->Fill(TMath::Log10(Ekin),1);
162                     } else if (Particle == 22) {
163                         Hits_g[nch-1]->Fill(r,(float) 1);
164                         HitDensity_g[nch-1]->Fill(r,wgt);
165                         Mom_g[nch-1]->Fill(TMath::Log10(Ekin),1);
166                     } else {
167                         Hits[nch-1]->Fill(r,(float) 1);
168                         HitDensity[nch-1]->Fill(r,wgt);
169                         Mom[nch-1]->Fill(TMath::Log10(Ekin),1);
170                     }
171                 } // hit loop
172             } // if MUON
173         } // track loop
174     }
175     
176 //Create a canvas, set the view range, show histograms
177     Int_t k;
178     TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
179     c1->Divide(2,4);
180     for (k=0; k<7; k++) {
181         c1->cd(k+1);
182         HitDensity[2*k]->Draw();
183     }
184
185     TCanvas *c2 = new TCanvas("c2","Hit Densities (gamma)",400,10,600,700);
186     c2->Divide(2,4);
187     for (k=0; k<7; k++) {
188         c2->cd(k+1);
189         HitDensity_g[2*k]->Draw();
190     }
191     
192     TCanvas *c3 = new TCanvas("c3","Hit Densities (neutron)",400,10,600,700);
193     c3->Divide(2,4);
194     for (k=0; k<7; k++) {
195         c3->cd(k+1);
196         HitDensity_n[2*k]->Draw();
197     }
198
199     TCanvas *c4 = new TCanvas("c4","Energy (charged)",400,10,600,700);
200     c4->Divide(2,4);
201     for (k=0; k<7; k++) {
202         c4->cd(k+1);
203         Mom[2*k]->Draw();
204     }
205
206     TCanvas *c5 = new TCanvas("c5","Energy (gamma)",400,10,600,700);
207     c5->Divide(2,4);
208     for (k=0; k<7; k++) {
209         c5->cd(k+1);
210         Mom_g[2*k]->Draw();
211     }
212
213     TCanvas *c6 = new TCanvas("c6","Energy (gamma)",400,10,600,700);
214     c6->Divide(2,4);
215     for (k=0; k<7; k++) {
216         c6->cd(k+1);
217         Mom_n[2*k]->Draw();
218     }
219 }
220
221
222
223
224
225
226