New RICH code J.Barbosa, A.Morsch, D.DiBari
[u/mrichter/AliRoot.git] / RICH / RICHtest.C
1 void RICHtest (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 // Connect the Root Galice file containing Geometry, Kine and Hits
17
18     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
19
20     
21     if (!file) {
22         printf("\n Creating galice.root \n");
23         file = new TFile("galice.root");
24     } else {
25         printf("\n galice.root found in file list");
26     }
27     file->ls();
28     
29 // Get AliRun object from file or create it if not on file
30     if (!gAlice) {
31         gAlice = (AliRun*)(file->Get("gAlice"));
32         if (gAlice) printf("AliRun object found on file\n");
33         if (!gAlice) {
34             printf("\n create new gAlice object");
35             gAlice = new AliRun("gAlice","Alice test program");
36         }
37     }
38     
39 //  Create some histograms
40
41
42    TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution"
43                            ,100,0.,500.);
44    TH1F *pcharge1 = new TH1F("pcharge1","Pad      Charge distribution"
45                            ,100,0.,200.);
46    TH1F *xresid1  = new TH1F("xresid1","x-Residuals"
47                            ,100,-0.5,0.5);
48    TH1F *yresid1  = new TH1F("yresid1","y-Residuals"
49                            ,100,-0.2,0.2);
50    TH1F *npads1   = new TH1F("npads1" ,"Pads per Hit"
51                            ,20,-0.5,19.5);
52    TH2F *xresys1  = new TH2F("xresys1","x-Residuals systematics"
53                            ,50,-2.0,2.0,100,-1.0,1.0);
54    TH2F *yresys1  = new TH2F("yresys1","y-Residuals systematics"
55                            ,50,-0.4,0.4,100,-0.1,0.1);
56
57    TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution"
58                            ,100,0.,500.);
59    TH1F *pcharge2 = new TH1F("pcharge2","Pad      Charge distribution"
60                            ,100,0.,200.);
61    TH1F *xresid2  = new TH1F("xresid2","x-Residuals"
62                            ,100,-0.5,0.5);
63    TH1F *yresid2  = new TH1F("yresid2","y-Residuals"
64                            ,100,-0.2,0.2);
65    TH1F *npads2   = new TH1F("npads2" ,"Pads per Hit"
66                            ,20,-0.5,19.5);
67    TH2F *xresys2  = new TH2F("xresys2","x-Residuals systematics"
68                            ,50,-2.0,2.0,100,-1.0,1.0);
69    TH2F *yresys2  = new TH2F("yresys2","y-Residuals systematics"
70                            ,50,-0.4,0.4,100,-0.1,0.1);
71
72
73    AliRICHchamber*  iChamber;
74 //
75 //   Loop over events 
76 //
77    Int_t Nh=0;
78    Int_t Nh1=0;
79    for (Int_t nev=0; nev<= evNumber2; nev++) {
80        //cout << "nev         " << nev <<endl;
81        printf ("Event Number       : %d\n",nev);
82        Int_t nparticles = gAlice->GetEvent(nev);
83        //cout << "nparticles  " << nparticles <<endl;
84        printf ("Number of particles: %d\n",nparticles);
85        if (nev < evNumber1) continue;
86        if (nparticles <= 0) return;
87
88        AliRICH *RICH  = (AliRICH*) gAlice->GetDetector("RICH");
89            printf("\n track %d %d \n ", nev, RICH);
90
91        TTree *TH = gAlice->TreeH();
92        Int_t ntracks = TH->GetEntries();
93        Int_t Nc=0;
94        Int_t   npad[2];
95        Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2];
96
97        printf ("Just got it");
98 //
99 //   Loop over events
100 //
101        for (Int_t track=0; track<ntracks;track++) {
102            
103            gAlice->ResetHits();
104            Int_t nbytes += TH->GetEvent(track);
105            if (RICH)  {
106                for(AliRICHhit* mHit=(AliRICHhit*)RICH->FirstHit(-1); 
107                    mHit;
108                    mHit=(AliRICHhit*)RICH->NextHit()) 
109                {
110                    Int_t   nch   = mHit->fChamber;  // chamber number
111                    Float_t x     = mHit->fX;        // x-pos of hit
112                    Float_t y     = mHit->fY;        // y-pos
113 //
114 //
115                    iChamber = & RICH->Chamber(nch-1);
116                    response=iChamber->GetResponseModel(mip);
117 //
118 //
119                    if (nch <= 1) {
120                        for (Int_t i=0; i<2; i++) {
121                            xmean[i]=0;
122                            ymean[i]=0;
123                            Q[i]=0;
124                            npad[i]=0;
125                        }
126                        for (AliRICHcluster *mPad=(AliRICHcluster*)RICH
127                                 ->FirstPad(mHit);
128                             mPad;
129                             mPad=(AliRICHcluster*)RICH->NextPad()
130                            )
131                        {
132                            Int_t nseg     = mPad->fCathode;   // segmentation 
133                            Int_t nhit     = mPad->fHitNumber; // hit number
134                            Int_t qtot     = mPad->fQ;         // charge
135                            Int_t ipx      = mPad->fPadX;      // pad number on X
136                            Int_t ipy      = mPad->fPadY;      // pad number on Y
137                            Int_t iqpad    = mPad->fQpad;      // charge per pad
138                            Int_t secpad   = mPad->fRSec;      // r-pos of pad
139 //
140 //
141                            segmentation=iChamber->GetSegmentationModel(nseg);
142                            Int_t ind=nseg-1;
143 //                     printf("Pad of hit, padx, pady, iqpad, ncha %d %d %d %d %d %d\n",
144 //                            nhit, ipx, ipy, iqpad, nseg, ind);
145
146                            
147                            if (nseg == 1) {
148                                pcharge1->Fill(iqpad,(float) 1);
149                                ccharge1->Fill(qtot ,(float) 1);
150                            }  else {
151                                pcharge2->Fill(iqpad,(float) 1);
152                                ccharge2->Fill(qtot ,(float) 1);
153                            }
154 // Calculate centre of gravity
155 //
156                            if (iqpad == 0 && ind==0) {
157                                printf("\n attention iqpad 0");
158                            }
159                            
160                            if (iqpad >  0 && secpad==1) {
161                                Float_t xc,yc;
162                                npad[ind]++;
163                                segmentation->GetPadCxy(ipx,ipy,xc,yc);
164 //                             printf("\n pad %d %d   ", ipx, ipy);
165 //                             printf("\n pos %f %f   ",  xc,  yc);
166                                xmean[ind]+=Float_t(iqpad*xc);
167                                ymean[ind]+=Float_t(iqpad*yc);
168                                Q[ind]+=iqpad;
169                            }
170                            
171                        } //pad hit loop
172                        for (Int_t i=0; i<2; i++) {
173                            segmentation = iChamber->GetSegmentationModel(i+1);
174                            
175 //                         if (npad[i] ==0) {
176 //                             printf("\n #pads=0 %f %f",x,y);
177 //                         }
178                            
179                            if (Q[i] >0) {
180                                xmean[i] = xmean[i]/Q[i];
181                                xres[i]  = xmean[i]-x;
182                                ymean[i] = ymean[i]/Q[i];
183                                yres[i]  = ymean[i]-y;
184 //                             printf("\n xmean, ymean %f %f",xmean[i],ymean[i]);
185 //                             printf("\n xhit, yhit %f %f",x,y);
186                                
187 // Systematic Error
188 //
189                                Int_t icx, icy;
190                                segmentation->
191                                    GetPadIxy(x,y,icx,icy);
192                                segmentation->
193                                    GetPadCxy(icx,icy,xonpad[i],yonpad[i]);
194                                xonpad[i]-=x;
195                                yonpad[i]-=y;
196                            }
197                        } // plane loop
198                        xresid1->Fill(xres[0],(float) 1);
199                        yresid1->Fill(yres[0],(float) 1);
200                        npads1->Fill(npad[0],(float) 1);
201                        if (npad[0] >=2 && Q[0] > 6 ) {
202                            xresys1->Fill(xonpad[0],xres[0],(float) 1);
203                            yresys1->Fill(yonpad[0],yres[0],(float) 1);
204                        }
205                        
206                        xresid2->Fill(xres[1],(float) 1);
207                        yresid2->Fill(yres[1],(float) 1);
208                        npads2->Fill(npad[1],(float) 1);
209                        if (npad[1] >=2 && Q[1] > 6) {
210                            xresys2->Fill(xonpad[1],xres[1],(float) 1);
211                            yresys2->Fill(yonpad[1],yres[1],(float) 1);
212                        }
213                    } // chamber 1
214                } // hit loop
215            } // if RICH
216        } // track loop
217    } // event loop 
218    
219 //Create a canvas, set the view range, show histograms
220    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
221    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
222    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
223    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
224    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
225    pad11->SetFillColor(11);
226    pad12->SetFillColor(11);
227    pad13->SetFillColor(11);
228    pad14->SetFillColor(11);
229    pad11->Draw();
230    pad12->Draw();
231    pad13->Draw();
232    pad14->Draw();
233    
234    pad11->cd();
235    ccharge1->SetFillColor(42);
236    ccharge1->SetXTitle("ADC units");
237    ccharge1->Draw();
238
239    pad12->cd();
240    pcharge1->SetFillColor(42);
241    pcharge1->SetXTitle("ADC units");
242    pcharge1->Draw();
243
244    pad13->cd();
245    xresid1->SetFillColor(42);
246    xresid1->Draw();
247
248    pad14->cd();
249    yresid1->SetFillColor(42);
250    yresid1->Draw();
251
252    TCanvas *c2 = new TCanvas("c2","Cluster Size",400,10,600,700);
253    pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
254    pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
255    pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
256    pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
257    pad21->SetFillColor(11);
258    pad22->SetFillColor(11);
259    pad23->SetFillColor(11);
260    pad24->SetFillColor(11);
261    pad21->Draw();
262    pad22->Draw();
263    pad23->Draw();
264    pad24->Draw();
265    
266    pad21->cd();
267    npads1->SetFillColor(42);
268    npads1->SetXTitle("Cluster Size");
269    npads1->Draw();
270
271    pad23->cd();
272    xresys1->SetXTitle("x on pad");
273    xresys1->SetYTitle("x-xcog");
274    xresys1->Draw();
275
276    pad24->cd();
277    yresys1->SetXTitle("y on pad");
278    yresys1->SetYTitle("y-ycog");
279    yresys1->Draw();
280
281    TCanvas *c3 = new TCanvas("c3","Charge and Residuals",400,10,600,700);
282    pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
283    pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
284    pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
285    pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
286    pad31->SetFillColor(11);
287    pad32->SetFillColor(11);
288    pad33->SetFillColor(11);
289    pad34->SetFillColor(11);
290    pad31->Draw();
291    pad32->Draw();
292    pad33->Draw();
293    pad34->Draw();
294    
295    pad31->cd();
296    ccharge2->SetFillColor(42);
297    ccharge2->SetXTitle("ADC units");
298    ccharge2->Draw();
299
300    pad32->cd();
301    pcharge2->SetFillColor(42);
302    pcharge2->SetXTitle("ADC units");
303    pcharge2->Draw();
304
305    pad33->cd();
306    xresid2->SetFillColor(42);
307    xresid2->Draw();
308
309    pad34->cd();
310    yresid2->SetFillColor(42);
311    yresid2->Draw();
312
313    TCanvas *c4 = new TCanvas("c4","Cluster Size",400,10,600,700);
314    pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
315    pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
316    pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
317    pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
318    pad41->SetFillColor(11);
319    pad42->SetFillColor(11);
320    pad43->SetFillColor(11);
321    pad44->SetFillColor(11);
322    pad41->Draw();
323    pad42->Draw();
324    pad43->Draw();
325    pad44->Draw();
326    
327    pad41->cd();
328    npads2->SetFillColor(42);
329    npads2->SetXTitle("Cluster Size");
330    npads2->Draw();
331
332    pad43->cd();
333    xresys2->SetXTitle("x on pad");
334    xresys2->SetYTitle("x-xcog");
335    xresys2->Draw();
336
337    pad44->cd();
338    yresys2->SetXTitle("y on pad");
339    yresys2->SetYTitle("y-ycog");
340    yresys2->Draw();
341    
342 }