]> git.uio.no Git - u/mrichter/AliRoot.git/blob - RICH/menu.C
Additional protection in case of zero magnetic field (Yu.Belikov)
[u/mrichter/AliRoot.git] / RICH / menu.C
1 AliRICH * R()    {return r;}
2 void ph(Int_t event=0) {R()->PrintHits(event);}    //utility print hits for 'event' event
3 void ps(Int_t event=0) {R()->PrintSDigits(event);} //utility print sdigits
4 void pd(Int_t event=0) {R()->PrintDigits(event);}  //utility print digits
5 void pc(Int_t event=0) {R()->PrintClusters(event);}//utility print clusters
6 void pt(Int_t event=0) {R()->PrintTracks(event);}  //utility print tracks
7
8 //__________________________________________________________________________________________________
9 void pp(int tid)
10 {
11   R()->GetLoader()->GetRunLoader()->LoadHeader();  R()->GetLoader()->GetRunLoader()->LoadKinematics();
12   
13   if(tid<0||tid>=R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack())
14     cout<<"Valid tid number is 0-"<<R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack()-1<<" for this event.\n";
15   else
16     PrintParticleInfo(tid);
17   
18   R()->GetLoader()->GetRunLoader()->UnloadKinematics();  R()->GetLoader()->GetRunLoader()->UnloadHeader();
19 }
20 //__________________________________________________________________________________________________
21 void PrintParticleInfo(int tid)
22 {
23 // Prints particle info for a given TID
24   TParticle *p=al->Stack()->Particle(tid);
25   cout<<p->GetName()<<"("<<tid<<")";
26   if(p->GetMother(0)!=-1){cout<<" from "; PrintParticleInfo(p->GetMother(0));}
27   else                   {cout<<endl;} 
28 }    
29 //__________________________________________________________________________________________________
30 Int_t prim(Int_t tid)
31 {
32 // Provides mother TID for the given TID
33   R()->GetLoader()->GetRunLoader()->LoadHeader();  R()->GetLoader()->GetRunLoader()->LoadKinematics();
34   
35   if(tid<0||tid>=R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack())
36     cout<<"Valid tid number is 0-"<<R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack()-1<<" for this event.\n";
37   else
38     while(1){
39       TParticle *p=R()->GetLoader()->GetRunLoader()->GetAliRun()->Stack()->Particle(tid);
40       if(p->GetMother(0)==-1) break;
41       tid=p->GetMother(0);
42     }
43   
44   R()->GetLoader()->GetRunLoader()->UnloadKinematics();  R()->GetLoader()->GetRunLoader()->UnloadHeader();
45   return tid;
46 }
47
48
49 //__________________________________________________________________________________________________
50
51 Double_t r2d = TMath::RadToDeg();
52 Double_t d2r = TMath::DegToRad();
53
54 void DisplFast(){ AliRICHDisplFast *d = new AliRICHDisplFast();  d->Exec();}  
55
56
57 void C_R()
58 {
59   AliRICHRecon *detect = new AliRICHRecon("RICH patrec algorithm","Reconstruction");
60     
61
62   for (int nev=0; nev< a->GetEventsPerRun(); nev++) {    // Event Loop
63     al->GetEvent(nev);
64     cout <<endl<< "Processing event:" <<nev<<endl;
65     detect->StartProcessEvent();
66   } // event loop  
67   delete detect;
68 }  
69 //__________________________________________________________________________________________________
70 void D_C()
71 {
72   TStopwatch sw;TDatime time;
73
74    AliRICHClusterFinder *z=new AliRICHClusterFinder(r); z->Exec();
75
76    cout << endl;
77    cout << "Info in Digits->Clusters: Start time: ";time.Print();
78    cout << "Info in Digits->Clusters: Stop  time: ";time.Set();  time.Print();
79    cout << "Info in Digits->Clusters: Time  used: ";sw.Print();
80 }
81 //__________________________________________________________________________________________________
82
83 void Show()
84 {  
85   cout<<endl;
86 //load all trees  
87   al->LoadHeader(); 
88     al->LoadKinematics();  
89       rl->LoadHits();  
90         Bool_t isSdigits=!rl->LoadSDigits();  
91           Bool_t isClusters=!rl->LoadRecPoints();
92             Bool_t isDigits=!rl->LoadDigits();//loaders
93   cout<<endl;  cout<<endl;  
94   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
95     Int_t iNparticles=a->GetEvent(iEventN);
96     Int_t iNprims=al->Stack()->GetNprimary();
97     
98     Int_t iTotalHits=0;
99     for(Int_t iPrimN=0;iPrimN<iNprims;iPrimN++){//prims loop
100       rl->TreeH()->GetEntry(iPrimN);      
101       iTotalHits+=r->Hits()->GetEntries();
102       TParticle *pPrim=al->Stack()->Particle(iPrimN);
103       if(iPrimN<10) Info("Show","Evt %4i prim %4i has %4i hits from %s (,%7.2f,%7.2f)",
104                   iEventN,iPrimN, r->Hits()->GetEntries(), pPrim->GetName(), pPrim->Theta()*r2d,pPrim->Phi()*r2d);
105     }//prims loop
106     Info("Show-HIT","Evt %i total:  %i particles %i primaries %i hits",
107                         iEventN,   iNparticles, iNprims,     iTotalHits);
108     if(isSdigits){
109       rl->TreeS()->GetEntry(0);
110       Info("Show-SDI","Evt %i contains %5i sdigits",iEventN,r->SDigits()->GetEntries());
111     }
112     if(isDigits){
113       rl->TreeD()->GetEntry(0);
114       for(int i=1;i<=7;i++)
115         Info("Show-DIG","Evt %i chamber %i contains %5i digits",
116                                  iEventN,   i,           r->Digits(i)->GetEntries());
117     }else
118         Info("Show-DIG","There is no digits for this event");
119     if(isClusters){
120       rl->TreeR()->GetEntry(0);
121       for(int i=1;i<=7;i++)
122         Info("Show-CLU","Evt %i chamber %i contains %5i clusters",
123                                  iEventN,   i,           r->Clusters(i)->GetEntries());
124     }
125     cout<<endl;
126   }//events loop
127 //unload all trees    
128   rl->UnloadHits();  
129     if(isSdigits) rl->UnloadSDigits(); 
130       if(isDigits) rl->UnloadDigits(); 
131         if(isClusters) rl->UnloadRecPoints();
132           al->UnloadHeader();
133             al->UnloadKinematics();
134   
135   TVector counters=r->Counters();
136   
137   counters(9)=counters(0); for(Int_t i=1;i<=8;i++) counters(9)-=counters(i);
138   counters.Print();
139   cout<<endl;
140 }//void Show()
141 //__________________________________________________________________________________________________
142
143
144 AliRun *a;
145 AliRunLoader *al;
146 AliLoader *rl,*tl,*il;
147
148 AliRICH *r;
149
150 Bool_t ReadAlice()
151 {
152   Info("ReadAlice","Tring to read ALICE from SIMULATED FILE.");
153   if(gAlice) delete gAlice;      
154   if(!(al=AliRunLoader::Open("galice.root","AlicE","update"))){
155     gSystem->Exec("rm -rf *.root *.dat");
156     Error("menu.C::ReadAlice","galice.root broken, removing all this garbage then init new one");
157     new AliRun("gAlice","Alice experiment system");
158     gAlice->Init("Config.C");
159     r=(AliRICH*)gAlice->GetDetector("RICH");
160     return kFALSE;
161   }
162   al->LoadgAlice();
163   if(!gAlice) Fatal("menu.C::ReadAlice","No gAlice in file");
164   a=al->GetAliRun();
165   a->SetDebug(0);    
166 //RICH      
167   if(!(r=(AliRICH*)gAlice->GetDetector("RICH"))) Warning("RICH/menu.C::ReadAlice","No RICH in file");
168   if(!(rl=al->GetLoader("RICHLoader")))          Warning("RICH/menu.C::ReadAlice","No RICH loader in file");        
169         
170   Info("ReadAlice","Run contains %i event(s)",gAlice->GetEventsPerRun());      
171   return kTRUE;
172 }
173 //__________________________________________________________________________________________________
174 void TestResponse()
175 {
176   TCanvas *pC=new TCanvas("c","Amplification test",900,800);
177   pC->Divide(1,2);
178   
179   
180   const Int_t nPoints=8;
181   THStack *pStackPhot=new THStack("StackPhot","photons");
182   THStack *pStackMip =new THStack("StackMip","mips");
183   TLegend *pLeg=new TLegend(0.6,0.2,0.9,0.5,"legend");    
184   TH1F *apHphot[nPoints];
185   TH1F *apHmip[nPoints];
186   
187   Double_t starty=0;
188   Double_t deltay=AliRICHParam::SectorSizeY()/nPoints;
189   
190   for(int i=0;i<nPoints;i++){
191     apHphot[i]=new TH1F(Form("hphot%i",i),"Qdc for Photon;QDC;Counts",500,0,500); apHphot[i]->SetLineColor(i);pStackPhot->Add(apHphot[i]);                 
192     apHmip[i] =new TH1F(Form("hmip%i",i),"Qdc for Mip;QDC;Counts",4000,0,4000);   apHmip[i]->SetLineColor(i);pStackMip->Add(apHmip[i]);                 
193     
194     pLeg->AddEntry(apHphot[i],Form("@(10,%5.2f->%5.2f)",starty+i*deltay,starty+i*deltay-AliRICHParam::SectorSizeY()/2));
195   }
196         
197   
198   TVector2 x2(0,0);  
199 //  AliRICHParam::ResetWireSag();
200   for(Int_t i=0;i<10000;i++){//events loop
201     for(int j=0;j<nPoints;j++){
202       x2.Set(10,starty+j*deltay);
203       apHphot[j]->Fill(AliRICHParam::TotQdc(x2,0));
204       apHmip[j]->Fill(AliRICHParam::TotQdc(x2,gRandom->Landau(600,150)*1e-9));
205     }
206   }
207   
208   pC->cd(1);  pStackMip->Draw("nostack");
209   pC->cd(2);  pStackPhot->Draw("nostack"); pLeg->Draw();
210 }//TestResponse()
211 //__________________________________________________________________________________________________
212 void TestSD()
213 {
214   Info("TestSD","Creating test sdigits.");
215   TVector3 hit(426.55,246.28,17.21);        
216   TVector2 x2=r->C(4)->Glob2Loc(hit);        
217   Int_t iTotQdc=r->P()->TotQdc(x2,624e-9);        
218   Int_t iPadXmin,iPadXmax,iPadYmin,iPadYmax;
219   Int_t padx,pady;
220   Int_t sec=r->P()->Loc2Pad(x2,padx,pady);
221   r->P()->Loc2Area(x2,iPadXmin,iPadYmin,iPadXmax,iPadYmax);
222   Info("TestSD","Initial hit (%7.2f,%7.2f,%7.2f)->(%7.2f,%7.2f)->(%4i,%4i,%4i) gives %i charge",
223                           hit.X(),hit.Y(),hit.Z(),x2.X(),x2.Y(),sec,padx,pady,iTotQdc);
224   
225   cout<<"left-down=("<<iPadXmin<<","<<iPadYmin<<") right-up=("<<iPadXmax<<','<<iPadYmax<<')'<<endl;
226   for(Int_t iPadY=iPadYmin;iPadY<=iPadYmax;iPadY++)
227     for(Int_t iPadX=iPadXmin;iPadX<=iPadXmax;iPadX++)
228        cout<<r->P()->FracQdc(x2,iPadX,iPadY)<<endl;
229   Info("TestSD","Stop.");
230 }//void TestSdigits()
231 //__________________________________________________________________________________________________
232 void TestC()
233 {
234   Info("TestC","Creating test clusters.");
235   R()->GetLoader()->MakeTree("R");R()->MakeBranch("R");
236   
237   AliRICHcluster c;
238   c.AddDigit(new AliRICHdigit(1,20,21,200,1,2,3));
239   c.AddDigit(new AliRICHdigit(1,22,21,250,1,2,3));
240   c.CoG();
241   
242   R()->AddCluster(c);  
243   
244   R()->GetLoader()->TreeR()->Fill();  R()->GetLoader()->WriteRecPoints("OVERWRITE");
245   R()->GetLoader()->UnloadRecPoints();
246   R()->ResetClusters();
247   
248   Info("TestC","Stop.");
249 }//TestC()
250 //__________________________________________________________________________________________________
251 void TestSeg()
252 {
253   AliRICHParam p;
254   
255   new TCanvas("name","PC segmentation");
256   gPad->Range(-20,-20,200,150);
257   AliRICHDisplFast::DrawSectors();
258   
259   TLatex t; t.SetTextSize(0.02);
260   t.DrawText(0,140,"View from interaction point");
261   t.DrawLatex(p.PcSizeX()+10,120,Form("Pc  %6.2fx%6.2fcm %3ix%3ipads",p.PcSizeX()    ,p.PcSizeY(),    p.NpadsX()   ,p.NpadsY()));
262   t.DrawLatex(p.PcSizeX()+10,115,Form("Sec %6.2fx%5.2fcm %3ix%2ipads",p.SectorSizeX(),p.SectorSizeY(),p.NpadsXsec(),p.NpadsYsec()));
263   t.DrawLatex(p.PcSizeX()+10,110,Form("Pad %6.2fx%4.2fcm DeadZone %6.2fcm",p.PadSizeX()   ,p.PadSizeY(),p.DeadZone()));
264   
265   TVector2 v2;
266   t.SetTextAlign(12);
267   v2=AliRICHParam::Pad2Loc( 1,24);  t.DrawText(v2.X(),v2.Y(),"sec 1");  
268   v2=AliRICHParam::Pad2Loc(81,24);  t.DrawText(v2.X(),v2.Y(),"sec 2");  
269   v2=AliRICHParam::Pad2Loc( 1,70);  t.DrawText(v2.X(),v2.Y(),"sec 3");  
270   v2=AliRICHParam::Pad2Loc(81,70);  t.DrawText(v2.X(),v2.Y(),"sec 4");  
271   v2=AliRICHParam::Pad2Loc( 1,120); t.DrawText(v2.X(),v2.Y(),"sec 5");  
272   v2=AliRICHParam::Pad2Loc(81,120); t.DrawText(v2.X(),v2.Y(),"sec 6");  
273   
274 //  TGaxis *pAx=new TGaxis(0,0,140,  0,0,140,510,"-="); pAx->SetTitle("x, cm"); pAx->SetTextSize(0.05); pAx->Draw();
275 //  TGaxis *pAy=new TGaxis(0,0,  0,140,0,140,510,"-="); pAy->SetTitle("y, cm"); pAy->SetTextSize(0.05); pAy->Draw();
276    
277   
278   t.SetTextColor(kBlue);  
279
280   Int_t padx,pady,sec;
281   Double_t margin=5;
282   Double_t x0=0; Double_t x1=p.SectorSizeX(); Double_t x2=p.SectorSizeX()+p.DeadZone(); Double_t x3=p.PcSizeX();
283   Double_t y0=0; Double_t y1=p.SectorSizeY(); Double_t y2=p.SectorSizeY()+p.DeadZone(); 
284   Double_t y3=2*p.SectorSizeY()+p.DeadZone(); Double_t y4=p.PcSizeY()-p.SectorSizeY();
285   Double_t y5=p.PcSizeY();
286    
287 //pads along x  
288   t.SetTextAlign(11); t.DrawText(x0,y5+margin,"1");  
289   t.SetTextAlign(31); t.DrawText(x1,y5+margin,"80");
290   t.SetTextAlign(11); t.DrawText(x2,y5+margin,"81");   
291   t.SetTextAlign(31); t.DrawText(x3,y5+margin,"160");
292 //pads along y    
293   t.SetTextAlign(11); t.DrawText(x3+margin,y0,"1");  
294   t.SetTextAlign(13); t.DrawText(x3+margin,y1,"48"); 
295   t.SetTextAlign(11); t.DrawText(x3+margin,y2,"49");
296   t.SetTextAlign(13); t.DrawText(x3+margin,y3,"96"); 
297   t.SetTextAlign(11); t.DrawText(x3+margin,y4,"97");
298   t.SetTextAlign(13); t.DrawText(x3+margin,y5,"144");        
299   
300   t.SetTextColor(kRed);
301    
302 //positions along x   
303   t.SetTextAlign(13);t.DrawText(x0,y0-margin,Form("%5.2f",x0));   
304   t.SetTextAlign(33);t.DrawText(x1,y0-margin,Form("%5.2f",x1));   
305   t.SetTextAlign(13);t.DrawText(x2,y0-margin,Form("%5.2f",x2));   
306   t.SetTextAlign(33);t.DrawText(x3,y0-margin,Form("%5.2f",x3));   
307 //positions along y
308   t.SetTextAlign(31);t.DrawText(x0-margin,y0,Form("%5.2f",y0));   
309   t.SetTextAlign(33);t.DrawText(x0-margin,y1,Form("%5.2f",y1));   
310   t.SetTextAlign(31);t.DrawText(x0-margin,y2,Form("%5.2f",y2));   
311   t.SetTextAlign(33);t.DrawText(x0-margin,y3,Form("%5.2f",y3));   
312   t.SetTextAlign(31);t.DrawText(x0-margin,y4,Form("%5.2f",y4));   
313   t.SetTextAlign(33);t.DrawText(x0-margin,y5,Form("%5.2f",y5));   
314 //coners      
315   t.SetTextSize(0.02);
316   TVector pad(2);
317   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x0,y0)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
318   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x0,y1)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
319   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x1,y1)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
320   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x1,y0)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
321
322   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x0,y2)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
323   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x0,y3)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
324   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x1,y3)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
325   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x1,y2)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
326    
327   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x0,y4)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
328   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x0,y5)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
329   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x1,y5)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
330   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x1,y4)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
331
332   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x2,y4)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
333   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x2,y5)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
334   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x3,y5)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
335   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x3,y4)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
336
337   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x2,y2)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
338   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x2,y3)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
339   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x3,y3)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
340   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x3,y2)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
341
342   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x2,y0)));t.SetTextAlign(11);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
343   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x2,y1)));t.SetTextAlign(13);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
344   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x3,y1)));t.SetTextAlign(33);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
345   v2=p.Pad2Loc(p.Loc2Pad(TVector2(x3,y0)));t.SetTextAlign(31);t.DrawText(v2.X(),v2.Y(),Form("%5.2f,%5.2f",v2.X(),v2.Y()));
346 //new canvas
347   new TCanvas("trasform","Test LRS-MRS transform");
348   
349   TView *pView=new TView(1);
350   pView->SetRange(-600,-600,-600,600,600,600);
351 //axis  
352   Double_t X[6]={0,0,0,300,0,0};  Double_t Y[6]={0,0,0,0,300,0};  Double_t Z[6]={0,0,0,0,0,300};  
353   TPolyLine3D *pXaxis=new TPolyLine3D(2,X);pXaxis->SetLineColor(kRed);   pXaxis->Draw();
354   TPolyLine3D *pYaxis=new TPolyLine3D(2,Y);pYaxis->SetLineColor(kGreen); pYaxis->Draw();
355   TPolyLine3D *pZaxis=new TPolyLine3D(2,Z);pZaxis->SetLineColor(kBlue);  pZaxis->Draw();  
356 //Draw PC for all chambers by trasfering Pc plane using Pc2Mrs methode  
357   Int_t iNpointsX=50,iNpointsY=50;  
358   for(Int_t iChamberN=1;iChamberN<=7;iChamberN++){//chamber loop
359     TPolyMarker3D *pChamber=new TPolyMarker3D(iNpointsX*iNpointsY);
360     Int_t i=0;
361     for(Double_t x=0;x<p.PcSizeX();x+=p.PcSizeX()/iNpointsX)
362       for(Double_t y=0;y<p.PcSizeY();y+=p.PcSizeY()/iNpointsY){//step loop
363         TVector3 v3=p.C(iChamberN)->Pc2Mrs(TVector2(x,y));//from regular grid of local PC points to MRS presentation
364         pChamber->SetPoint(i++,v3.X(),v3.Y(),v3.Z());//Pc plane poing in MRS
365       }//step loop
366     pChamber->SetMarkerSize(1);
367     pChamber->SetMarkerColor(iChamberN);
368     pChamber->Draw();  
369     t.SetNDC();t.SetTextColor(iChamberN); t.DrawText(0.1,iChamberN*0.1,Form("Chamber %i",iChamberN));    
370   }//chamber loop   
371 //  gPad->GetView()->RotateView(94,45);
372 }//void TestSeg()
373 //__________________________________________________________________________________________________
374 void TestMenu()
375 {
376   TControlBar *pMenu = new TControlBar("vertical","RICH test");
377   pMenu->AddButton("Test segmentation",  "TestSeg()",         "Test AliRICHParam segmentation methods");
378   pMenu->AddButton("Test response",      "TestResponse()",    "Test AliRICHParam response methods");
379   pMenu->AddButton("Test sdigits",       "TestSD()",          "Create test set of sdigits");
380   pMenu->AddButton("Test clusters",      "TestC()",           "Create test set of clusters");
381   pMenu->Show();  
382 }//TestMenu()
383 //__________________________________________________________________________________________________
384 void menu()
385
386   TControlBar *pMenu = new TControlBar("vertical","RICH main");
387        
388   if(ReadAlice()){//it's from file, reconstruct
389     pMenu->AddButton("Show",            "Show()",             "Shows the structure of events in files");
390     pMenu->AddButton("Display Fast",    "DisplFast()",        "Display Fast");
391     pMenu->AddButton("Control Plots",   "R()->ControlPlots()","Create some control histograms");
392     
393   }else{//it's aliroot, simulate
394     pMenu->AddButton("Run",         "a->Run(1)",       "Process!");
395     pMenu->AddButton("Geo GUI", "new G3GeometryGUI;","Create instance of G4GeometryGUI"); 
396     pMenu->AddButton("Read RAW","ReadRaw()","Read a list of digits from test beam file"); 
397   }
398   pMenu->AddButton("Test submenu",    "TestMenu()",            "Shows test submenu");
399   pMenu->AddButton("Browser",         "new TBrowser;",         "Start ROOT TBrowser");
400   pMenu->AddButton("Debug ON",     "DebugON();",   "Switch debug on-off");   
401   pMenu->AddButton("Debug OFF",    "DebugOFF();",   "Switch debug on-off");   
402   pMenu->AddButton("Quit",            ".q",                    "Close session");
403   pMenu->Show();
404   a=gAlice;//for manual manipulation convinience
405 }//menu()
406 //__________________________________________________________________________________________________
407 void DebugOFF(){  Info("DebugOFF","");  a->SetDebug(0);  r->SetDebug(0);  AliLoader::SetDebug(0);}
408 void DebugON() {  Info("DebugON","");   a->SetDebug(1);  r->SetDebug(1);  AliLoader::SetDebug(1);}
409
410 void ReadRaw()
411 {
412 // First interation for TB raw data reader  
413   char *sRawFileName="beam/run1822_4sigma.dat";
414   FILE *pRawFile=fopen(sRawFileName,"r");
415   if(!pRawFile){Error("ReadRaw","Something wrong with raw data file %s",sRawFileName); return;}
416   
417   Int_t iNpads=0,q=0,x=0,y=0;
418   
419   while(fscanf(pRawFile,"%i",&iNpads)){
420     Info("ReadRaw","%i pads fired",iNpads);
421     for(Int_t i=0;i<iNpads;i++)
422       fscanf(pRawFile,"%i %i %i",&q,&x,&y);
423   }
424   fclose(pRawFile);
425 }
426 //__________________________________________________________________________________________________
427 void FillContribs(Int_t part,Int_t C, Bool_t print)
428 {
429   static Int_t contrib[10][7][2];
430   static Bool_t kEnter=kFALSE;
431
432   C--; // chamber numbering from 1 to 7
433   if(!kEnter) {
434     for(Int_t i=0;i<10;i++) {
435       for(Int_t j=0;j<7;j++) {
436         for(Int_t k=0;k<2;k++) contrib[i][j][k]=0;
437       }
438     }
439     kEnter=kTRUE;
440   }
441
442   if(print) {
443     for(Int_t k=0;k<2;k++) {cout << "----------------Positives  to RICH ---------------" << endl;
444                             cout << " chamber    1    2    3     4     5     6    7    " << endl;
445                             cout << " -------------------------------------------------" << endl;
446       for(Int_t i=0;i<4;i++) {
447         if(i==0) cout << " e";
448         if(i==1) cout << "pi";
449         if(i==2) cout << " K";
450         if(i==3) cout << " p";
451         if(k==0) cout << "+         ";
452         if(k==1) cout << "-         ";
453         for(Int_t j=0;j<7;j++) {
454           cout << contrib[i][j][k] << "    ";
455         }
456           cout << endl;
457       }
458     }
459   }
460
461   // +ves
462   if(part==kPositron)    contrib[0][C][0]++;
463   if(part==kPiPlus)      contrib[1][C][0]++;
464   if(part==kKPlus)       contrib[2][C][0]++;
465   if(part==kProton)      contrib[3][C][0]++;
466
467   // -ves
468   if(part==kElectron)    contrib[0][C][1]++;
469   if(part==kPiMinus)     contrib[1][C][1]++;
470   if(part==kKMinus)      contrib[2][C][1]++;
471   if(part==kProtonBar)   contrib[3][C][1]++;
472 }
473 //__________________________________________________________________________________________________
474 void ParticleContribs()
475 {
476
477   TH1F *distH1 = new TH1F("distH1","distH1",100,0.,5.);
478   AliRICH *pRich = (AliRICH*)gAlice->GetDetector("RICH");
479   Bool_t isHits    =!pRich->GetLoader()->LoadHits();
480
481   pRich->GetLoader()->GetRunLoader()->LoadHeader();  pRich->GetLoader()->GetRunLoader()->LoadKinematics();
482
483
484   if(!isHits){Error("Exec","No hits. No contribs to calculate.");return;}
485
486   for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events Loop
487     pRich->GetLoader()->GetRunLoader()->GetEvent(iEventN);
488     cout << " event n. " << iEventN << endl;
489     Int_t nPrimaries = (Int_t)pRich->GetLoader()->TreeH()->GetEntries();
490     TObjArray * Hits = new TObjArray[nPrimaries];
491
492     for(Int_t i=0;i<nPrimaries;i++) {
493       pRich->GetLoader()->TreeH()->GetEntry(i);
494       Int_t nHits = pRich->Hits()->GetEntries();
495       for(Int_t k=0;k<nHits;k++)         Hits[i].Add(pRich->Hits()->At(k));
496
497     }
498 //deals with hits
499       for(Int_t i=0;i<nPrimaries;i++){//prims loop
500         pRich->GetLoader()->TreeH()->GetEntry(i);
501         Int_t nHits = pRich->Hits()->GetEntries();
502         for(Int_t j=0;j<nHits;j++){//hits loop
503           AliRICHhit *pHit = (AliRICHhit*)Hits[i].At(j);
504           TParticle *pParticle = pRich->GetLoader()->GetRunLoader()->GetAliRun()->Stack()->Particle(pHit->GetTrack());
505 //          if(pParticle->P()>1) FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE);
506           FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE);
507           if(pParticle->GetPDG()->Charge()!=0) distH1->Fill(pHit->Length());
508         }//hits loop
509       }//prims loop
510     }// event loop
511     FillContribs(0,0,kTRUE);
512     distH1->Draw();
513
514 }//ParticleContribs()