get tables from the aliroot directory if they are not in the current one
[u/mrichter/AliRoot.git] / MUON / MUONhitrectest.C
CommitLineData
a9e2aefa 1void MUONhitrectest (Int_t evNumber1=0,Int_t evNumber2=0)
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
12 if (gClassTable->GetID("AliRun") < 0) {
13 gROOT->LoadMacro("loadlibs.C");
14 loadlibs();
15 }
16
17// Connect the Root Galice file containing Geometry, Kine and Hits
18
19 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
20 if (!file) file = new TFile("galice.root");
21
22// Get AliRun object from file or create it if not on file
23
24 if (!gAlice) {
25 gAlice = (AliRun*)file->Get("gAlice");
26 if (gAlice) printf("AliRun object found on file\n");
27 if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
28 }
29// Create some histograms
30
31 TH2F *h21 = new TH2F("h21","Hits",100,-100,100,100,-100,100);
32 TH2F *h22 = new TH2F("h22","CoG ",100,-100,100,100,-100,100);
33 TH1F *h1 = new TH1F("h1","Multiplicity",30,-0.5,29.5);
34 TH1F *hmult = new TH1F("hmult","Multiplicity",30,-0.5,29.5);
35 TH1F *hresx = new TH1F("hresx","Residuals",100,-1,1);
36 TH1F *hresy = new TH1F("hresy","Residuals",100,-10.,10);
37 TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
38 TH1F *hpos = new TH1F("hnpos","Possibilities",10,-0.5,9.5);
39 TH2F *hchi1 = new TH2F("hchi1","Chi2 vs Residuals",100,0,0.2,100,-500,500);
40 TH2F *hchi2 = new TH2F("hchi2","Chi2 vs Residuals",100,0,20,100,-500,500);
41//
42// Loop over events
43//
44 Int_t Nh=0;
45 Int_t Nh1=0;
46
47 for (int nev=0; nev<= evNumber2; nev++) {
48 Int_t nparticles = gAlice->GetEvent(nev);
49 cout << "nev " << nev <<endl;
50 cout << "nparticles " << nparticles <<endl;
51 if (nev < evNumber1) continue;
52 if (nparticles <= 0) return;
53
54 TTree *TH = gAlice->TreeH();
55 Int_t ntracks = TH->GetEntries();
56 cout<<"ntracks "<<ntracks<<endl;
57
58 Int_t nbytes = 0;
59
60
61
62// Get pointers to Alice detectors and Digits containers
63 AliMUON *MUON = (AliMUON*) gAlice->GetModule("MUON");
64 TClonesArray *Particles = gAlice->Particles();
65 TTree *TR = gAlice->TreeR();
66 Int_t nent=TR->GetEntries();
4988e518 67
a9e2aefa 68 if (MUON) {
69 Float_t xc[2000], yc[2000], itrc[2000];
70 TClonesArray *MUONrawclust = MUON->RawClustAddress(0);
4988e518 71 nbytes += TR->GetEvent(0);
a9e2aefa 72 Int_t nrawcl = MUONrawclust->GetEntries();
73 AliMUONRawCluster *mRaw;
74 for (Int_t iraw=0; iraw < nrawcl; iraw++) {
75 mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
76 xc[iraw]=mRaw->fX[1];
77 yc[iraw]=mRaw->fY[0];
78 itrc[iraw]=mRaw->fTracks[1];
79 } // cluster
80
81 gAlice->ResetHits();
82 for (Int_t itrack=0; itrack<ntracks ;itrack++)
83 {
84 Int_t nbytes += TH->GetEvent(itrack);
85 Int_t nhit=MUON->Hits()->GetEntriesFast();
86 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(itrack);
87 mHit;
88 mHit=(AliMUONHit*)MUON->NextHit())
89 {
4988e518 90 Int_t nch = mHit->Chamber(); // chamber number
91 Float_t x = mHit->X(); // x-pos of hit
92 Float_t y = mHit->Y(); // y-pos
93 Int_t ip = mHit->Particle();
a9e2aefa 94 if (ip != kMuonPlus && ip != kMuonMinus) continue;
95 if (nch != 1) continue;
96 Float_t dmin=1.e8;
97 Int_t imin=-1;
98//
99// Loop over clusters
100 Int_t npos=0;
101
102 for (Int_t i=0; i<nrawcl; i++) {
103 Float_t dx = xc[i]-x;
104 Float_t dy = yc[i]-y;
4988e518 105 Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
a9e2aefa 106 if (dr<dmin) {
107 imin=i;
108 dmin=dr;
109 }
4988e518 110 if (TMath::Abs(dy) < 0.05 && TMath::Abs(dx) < 0.2 ) npos++;
a9e2aefa 111 }
112 mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(imin);
4988e518 113 Int_t track = mRaw->fTracks[1];
114 Float_t xrec = mRaw->fX[1];
115 Float_t yrec = mRaw->fY[0];
116 Float_t x_res = xrec-x;
117 Float_t y_res = yrec-y;
118 if (y_res > 0.05 ) printf("\n track %d %f\n", itrack, y_res);
119
a9e2aefa 120 hresx->Fill(x_res,1.);
121 hresy->Fill(y_res,1.);
122 hpos->Fill(Float_t(npos), 1.);
123 hresym->Fill(y_res*1.e4,1.);
124 if (npos == 0) {
125 printf("No cluster %d %f %f %d\n",
126 itrack, x, y, nhit);
127 h21->Fill(x,y,1.);
128 }
a9e2aefa 129 } // hits
130 } // tracks
131 } // end if MUON
132 } // event loop
133 TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
134 c1->Divide(2,2);
135 c1->cd(1);
136 hresx->SetFillColor(42);
137 hresx->SetXTitle("xrec-x");
138 hresx->Draw();
139
140 c1->cd(2);
141 hresy->SetFillColor(42);
142 hresy->SetXTitle("yrec-y");
143 hresy->Draw();
144
145 c1->cd(3);
146 hpos->SetFillColor(42);
147 hpos->SetXTitle("Possibilities");
148 hpos->Draw();
149
150 c1->cd(4);
151 hresym->SetFillColor(42);
152 hresym->SetXTitle("yrec-y");
153 hresym->Fit("gaus");
154 hresym->Draw();
155
156 TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
157 c2->Divide(2,2);
158
159 c2->cd(1);
160 h21->SetFillColor(42);
161 h21->SetXTitle("x");
162 h21->SetYTitle("y");
163 h21->Draw();
164
165 c2->cd(2);
166 hmult->SetFillColor(42);
167 hmult->SetXTitle("multiplicity");
168 hmult->Draw();
169
170 c2->cd(3);
171 hchi1->SetFillColor(42);
172 hchi1->SetXTitle("chi2");
173 hchi1->Draw();
174
175 c2->cd(4);
176 hchi2->SetFillColor(42);
177 hchi2->SetXTitle("chi2");
178 hchi2->Draw();
179
180}
181
182
183