Increment the version number
[u/mrichter/AliRoot.git] / RICH / RICHgainvar.C
1 void RICHgainvar (Int_t evNumber1=0,Int_t evNumber2=0) 
2 {
3
4   gClassTable->GetID("AliRun");
5
6
7   // Dynamically link some shared libs
8  
9   if (gClassTable->GetID("AliRun") < 0) {
10     gROOT->LoadMacro("loadlibs.C");
11     loadlibs();
12   }else {
13     delete gAlice;
14     gAlice = 0;
15   }
16   
17   gAlice=0;
18   
19   // Connect the Root Galice file containing Geometry, Kine and Hits
20   
21   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
22   if (!file) file = new TFile("galice.root","UPDATE");
23   
24   // Get AliRun object from file or create it if not on file
25   
26   if (!gAlice) {
27     gAlice = (AliRun*)file->Get("gAlice");
28     if (gAlice) printf("AliRun object found on file\n");
29     if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
30   }
31   else {
32     delete gAlice;
33     gAlice = (AliRun*)file->Get("gAlice");
34     if (gAlice) printf("AliRun object found on file\n");
35     if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
36   }
37   
38   AliRICH *RICH  = (AliRICH*)gAlice->GetDetector("RICH");
39   
40   TH1F *gainX = new TH1F("gainX","Gain along the wires",40,0,40);
41   TH1F *gainIn = new TH1F("gainIn","Gain near the center of photocathode",300,0,300);
42   TH1F *gainOut = new TH1F("gainOut","Gain near the edge of photocathode",300,0,300);
43   /*TH1F *gainPhi1 = new TH1F("gainPhi1","Gain variation along wires, 0-30 degrees",300,0,300);
44   TH1F *gainPhi2 = new TH1F("gainPhi2","Gain variation along wires, 30-60 degrees",300,0,300);
45   TH1F *gainPhi3 = new TH1F("gainPhi3","Gain variation along wires, 60-90 degrees",300,0,300);
46   TH1F *gainPhi4 = new TH1F("gainPhi4","Gain variation along wires, 90-120 degrees",300,0,300);
47   TH1F *gainPhi5 = new TH1F("gainPhi5","Gain variation along wires, 120-150 degrees",300,0,300);
48   TH1F *gainPhi6 = new TH1F("gainPhi6","Gain variation along wires, 150-180 degrees",300,0,300);
49   TH1F *gainPhi7 = new TH1F("gainPhi7","Gain variation along wires, 180-210 degrees",300,0,300);
50   TH1F *gainPhi8 = new TH1F("gainPhi8","Gain variation along wires, 210-240 degrees",300,0,300);
51   TH1F *gainPhi9 = new TH1F("gainPhi9","Gain variation along wires, 240-270 degrees",300,0,300);
52   TH1F *gainPhi10 = new TH1F("gainPhi10","Gain variation along wires, 270-300 degrees",300,0,300);
53   TH1F *gainPhi11 = new TH1F("gainPhi11","Gain variation along wires, 300-330 degrees",300,0,300);
54   TH1F *gainPhi12 = new TH1F("gainPhi12","Gain variation along wires, 330-360 degrees",300,0,300);
55   TH1F *gainPhi = new TH1F("gainPhi","Gain variation along phi",36,0,360);*/
56
57
58   TH1F *gainPhi1 = new TH1F("gainPhi1","Gain variation along wires, 315 to 45 degrees",300,0,300);
59   TH1F *gainPhi2 = new TH1F("gainPhi2","Gain variation along wires, 45 to 135 degrees",300,0,300);
60   TH1F *gainPhi3 = new TH1F("gainPhi3","Gain variation along wires, 135 to 225 degrees",300,0,300);
61   TH1F *gainPhi4 = new TH1F("gainPhi4","Gain variation along wires, 225 to 315 degrees",300,0,300);
62   TH1F *gainPhi = new TH1F("gainPhi","Gain variation along phi",8,0,360);
63   
64   //   Start loop over events 
65
66   for (int nev=0; nev<= evNumber2; nev++) {
67     Int_t nparticles = gAlice->GetEvent(nev);
68     
69     
70     printf ("\n**********************************\nProcessing Event: %d\n",nev);
71     printf ("Particles       : %d\n\n",nparticles);
72     if (nev < evNumber1) continue;
73     if (nparticles <= 0) return;
74        
75     // Get pointers to RICH detector and Hits containers
76        
77     
78     TTree *TH = gAlice->TreeH(); 
79     Int_t ntracks = TH->GetEntries();
80     Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
81     gAlice->TreeR()->GetEvent(nent-1);
82     TClonesArray *Rawclusters = RICH->RawClustAddress(2);    //  Raw clusters branch
83     //printf ("Rawclusters:%p",Rawclusters);
84     Int_t nrawclusters = Rawclusters->GetEntriesFast();
85     //printf (" nrawclusters:%d\n",nrawclusters);
86     
87     
88     // Start loop on tracks in the hits containers
89     Int_t Nc=0;
90     for (Int_t track=0; track<ntracks;track++) {
91       printf ("\nProcessing Track: %d\n",track);
92       gAlice->ResetHits();
93       Int_t nbytes += TH->GetEvent(track);
94       if (RICH)
95         TClonesArray *Hits = RICH->Hits();            // Hits branch
96       //see hits distribution
97       
98       Int_t nhits = Hits->GetEntriesFast();
99       printf("Hits            : %d\n",nhits);
100       for (Int_t hit=0;hit<nhits;hit++) {
101         mHit = (AliRICHHit*) Hits->UncheckedAt(hit);
102         Float_t mx  = mHit->X();                    // x-pos of hit
103         Float_t my  = mHit->Z();                    // y-pos
104       }
105
106
107       if (nrawclusters) {
108         printf("Raw Clusters    : %d\n",nrawclusters);
109         for (Int_t hit=0;hit<nrawclusters;hit++) {
110           rcHit = (AliRICHRawCluster*) Rawclusters->UncheckedAt(hit);
111           //Int_t nchamber = rcHit->fChamber;     // chamber number
112           //Int_t nhit = cHit->fHitNumber;        // hit number
113           Int_t qtot = rcHit->fQ;                 // charge
114           Int_t fx  =  rcHit->fX;                 // x-position
115           Int_t fy  =  rcHit->fY;                 // y-position
116           Int_t mult = rcHit->fMultiplicity;      // How many pads form the cluster
117                             
118           Float_t radfid = TMath::Sqrt((TMath::Power(fx-mx,2)+(TMath::Power(fy-my,2))));
119           //printf("Radius:%f\n", radfid);
120           
121           
122           if(qtot<200 && radfid>=9 && radfid<=15)
123             {
124               if (fy>=0)
125                 gainX->Fill(fy-40,(float) qtot);
126               else
127                 gainX->Fill(fy+40,(float) qtot);
128               
129               if(fy<=0 && fy>=-4)
130                 gainOut->Fill(qtot,(float) 1);
131               if(fy<=-11 && fy>=-15)
132                 gainIn->Fill(qtot,(float) 1);
133               
134               Float_t phi = TMath::ATan2(fx-mx, fy-my);
135               Float_t phi_deg = phi/TMath::Pi()*180;
136               
137               if ((fx-mx)<0)
138                 {
139                   phi_deg = 360+phi_deg;
140                   //printf("X:%f, Y:%f, Phi:%f\n",-fy+my,fx-mx,phi_deg);
141                 }
142               else
143                 {
144                   //printf("X:%f, Y:%f, Phi:%f\n",-fy+my,fx-mx,phi_deg);
145                 }
146               
147               
148               /*if (phi_deg>=0 && phi_deg<30)
149                 gainPhi1->Fill(qtot,(float) 1);
150               if (phi_deg>=30 && phi_deg<60)
151                 gainPhi2->Fill(qtot,(float) 1);
152               if (phi_deg>=60 && phi_deg<90)
153                 gainPhi3->Fill(qtot,(float) 1);
154               if (phi_deg>=90 && phi_deg<120)
155                 gainPhi4->Fill(qtot,(float) 1);
156               if (phi_deg>=120 && phi_deg<150)
157                 gainPhi5->Fill(qtot,(float) 1);
158               if (phi_deg>=150 && phi_deg<180)
159                 gainPhi6->Fill(qtot,(float) 1);
160               if (phi_deg>=180 && phi_deg<210)
161                 gainPhi7->Fill(qtot,(float) 1);
162               if (phi_deg>=210 && phi_deg<240)
163                 gainPhi8->Fill(qtot,(float) 1);
164               if (phi_deg>=240 && phi_deg<270)
165                 gainPhi9->Fill(qtot,(float) 1);
166               if (phi_deg>=270 && phi_deg<300)
167                 gainPhi10->Fill(qtot,(float) 1);
168               if (phi_deg>=300 && phi_deg<330)
169                 gainPhi11->Fill(qtot,(float) 1);
170               if (phi_deg>=330 && phi_deg<360)
171                 gainPhi12->Fill(qtot,(float) 1);*/
172               
173               if (phi_deg>=315 || phi_deg<45)
174                 gainPhi1->Fill(qtot,(float) 1);
175               if (phi_deg>=45 && phi_deg<135)
176                 gainPhi2->Fill(qtot,(float) 1);
177               if (phi_deg>=135 && phi_deg<225)
178                 gainPhi3->Fill(qtot,(float) 1);
179               if (phi_deg>=225 && phi_deg<315)
180                 gainPhi4->Fill(qtot,(float) 1);
181               
182             }
183           
184           //printf("Y:%d, Q:%d %d in %d\n",fy, qtot,hit,nrawclusters);
185         }
186       }
187     }
188   }
189
190   
191   TCanvas *c16 = new TCanvas("c16","Gain Variation over Phi",50,50,1200,700);
192   //c16->Divide(4,2);
193   //c16->SetFillColor(42);
194   
195   c16->Divide(2,2);
196   
197   c16->cd(1);
198   c16_1->SetLogy();
199   gainPhi1->SetFillColor(5);
200   gainPhi1->SetXTitle("ADC counts");
201   gainPhi1->Fit("expo","M");
202   expo->SetLineColor(2);
203   expo->SetLineWidth(3);
204   Float_t inv_slope1 = -1/expo->GetParameter(1);
205   //gainPhi->Fill(15, (float) inv_slope1);
206   gainPhi->Fill(0, (float) inv_slope1);
207   gainPhi1->Draw();
208   
209   c16->cd(2);
210   c16_2->SetLogy();
211   gainPhi2->SetFillColor(5);
212   gainPhi2->SetXTitle("ADC counts");
213   gainPhi2->Fit("expo","M");
214   expo->SetLineColor(2);
215   expo->SetLineWidth(3);
216   Float_t inv_slope2 = -1/expo->GetParameter(1);
217   //gainPhi->Fill(45, (float) inv_slope2);
218   gainPhi->Fill(90, (float) inv_slope2);
219   gainPhi2->Draw();
220   
221   c16->cd(3);
222   c16_3->SetLogy();
223   gainPhi3->SetFillColor(5);
224   gainPhi3->SetXTitle("ADC counts");
225   gainPhi3->Fit("expo","M");
226   expo->SetLineColor(2);
227   expo->SetLineWidth(3);
228   Float_t inv_slope3 = -1/expo->GetParameter(1);
229   //gainPhi->Fill(75, (float) inv_slope3);
230   gainPhi->Fill(180, (float) inv_slope3);
231   gainPhi3->Draw();
232   
233   c16->cd(4);
234   c16_4->SetLogy();
235   gainPhi4->SetFillColor(5);
236   gainPhi4->SetXTitle("ADC counts");
237   gainPhi4->Fit("expo","M");
238   expo->SetLineColor(2);
239   expo->SetLineWidth(3);
240   Float_t inv_slope4 = -1/expo->GetParameter(1);
241   //gainPhi->Fill(105, (float) inv_slope4);
242   gainPhi->Fill(270, (float) inv_slope4);
243   gainPhi4->Draw();
244   
245   /*c16->cd(5);
246   c16_5->SetLogy();
247   gainPhi5->SetFillColor(5);
248   gainPhi5->SetXTitle("ADC counts");
249   gainPhi5->Fit("expo","M");
250   expo->SetLineColor(2);
251   expo->SetLineWidth(3);
252   Float_t inv_slope5 = -1/expo->GetParameter(1);
253   gainPhi->Fill(135, (float) inv_slope5);
254   gainPhi5->Draw();
255   
256   c16->cd(6);
257   c16_6->SetLogy();
258   gainPhi6->SetFillColor(5);
259   gainPhi6->SetXTitle("ADC counts");
260   gainPhi6->Fit("expo","M");
261   expo->SetLineColor(2);
262   expo->SetLineWidth(3);
263   Float_t inv_slope6 = -1/expo->GetParameter(1);
264   gainPhi->Fill(165, (float) inv_slope6);
265   gainPhi6->Draw();
266   
267   c16->cd(7);
268   c16_7->SetLogy();
269   gainPhi7->SetFillColor(5);
270   gainPhi7->SetXTitle("ADC counts");
271   gainPhi7->Fit("expo","M");
272   expo->SetLineColor(2);
273   expo->SetLineWidth(3);
274   Float_t inv_slope7 = -1/expo->GetParameter(1);
275   gainPhi->Fill(195, (float) inv_slope7);
276   gainPhi7->Draw();
277   
278   c16->cd(8);
279   c16_8->SetLogy();
280   gainPhi8->SetFillColor(5);
281   gainPhi8->SetXTitle("ADC counts");
282   gainPhi8->Fit("expo","M");
283   expo->SetLineColor(2);
284   expo->SetLineWidth(3);
285   Float_t inv_slope8 = -1/expo->GetParameter(1);
286   gainPhi->Fill(225, (float) inv_slope8);
287   gainPhi8->Draw();
288   
289   c16->cd(9);
290   c16_9->SetLogy();
291   gainPhi9->SetFillColor(5);
292   gainPhi9->SetXTitle("ADC counts");
293   gainPhi9->Fit("expo","M");
294   expo->SetLineColor(2);
295   expo->SetLineWidth(3);
296   Float_t inv_slope9 = -1/expo->GetParameter(1);
297   gainPhi->Fill(255, (float) inv_slope9);
298   c16->Update();
299   gainPhi9->Draw();
300   
301   
302   c16->cd(10);
303   c16_10->SetLogy();
304   gainPhi10->SetFillColor(5);
305   gainPhi10->SetXTitle("ADC counts");
306   gainPhi10->Fit("expo","M");
307   expo->SetLineColor(2);
308   expo->SetLineWidth(3);
309   Float_t inv_slope10 = -1/expo->GetParameter(1);
310   gainPhi->Fill(285, (float) inv_slope10);
311   gainPhi10->Draw();
312   
313   c16->cd(11);
314   c16_11->SetLogy();
315   gainPhi11->SetFillColor(5);
316   gainPhi11->SetXTitle("ADC counts");
317   gainPhi11->Fit("expo","M");
318   expo->SetLineColor(2);
319   expo->SetLineWidth(3);
320   Float_t inv_slope11 = -1/expo->GetParameter(1);
321   gainPhi->Fill(315, (float) inv_slope11);
322   gainPhi11->Draw();
323   
324   c16->cd(12);
325   c16_12->SetLogy();
326   gainPhi12->SetFillColor(5);
327   gainPhi12->SetXTitle("ADC counts");
328   gainPhi12->Fit("expo","M");
329   expo->SetLineColor(2);
330   expo->SetLineWidth(3);
331   Float_t inv_slope12 = -1/expo->GetParameter(1);
332   gainPhi->Fill(335, (float) inv_slope12);
333   gainPhi12->Draw();*/
334   
335   TCanvas *c17 = new TCanvas("c17","Gain Variation",50,50,600,700);
336   c17->Divide(2,2);
337   
338   //c11->SetFillColor(42);
339   
340   c17->cd(1);
341   c17_1->SetLogy();
342   gainIn->SetFillColor(5);
343   gainIn->SetXTitle("ADC counts");
344   gainIn->Fit("expo","M");
345   expo->SetLineColor(2);
346   expo->SetLineWidth(3);
347   gainIn->Draw();
348   
349   c17->cd(2);
350   c17_2->SetLogy();
351   gainOut->SetFillColor(5);
352   gainOut->SetXTitle("ADC counts");
353   gainOut->Fit("expo","M");
354   expo->SetLineColor(2);
355   expo->SetLineWidth(3);
356   gainOut->Draw();
357   
358   c17->cd(3);
359   gainX->SetFillColor(5);
360   gainX->SetYTitle("Total ADC counts");
361   gainX->SetXTitle("(pads)");
362   gainX->Draw();
363   
364   c17->cd(4);
365   gainPhi->SetFillColor(5);
366   gainPhi->SetYTitle("Gain (ADC)");
367   gainPhi->SetXTitle("phi (degrees)");
368   gainPhi->Draw();
369
370   printf("\nEnd of macro\n");
371   printf("**********************************\n");
372   
373 }