4ace753ff23c1b40b3e390d44e8906e106bbb918
[u/mrichter/AliRoot.git] / RICH / AliRICHDetect.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 #include <stdlib.h>
19
20
21 #include "AliRICH.h"
22 #include "AliRICHPoints.h"
23 #include "AliRICHDetect.h"
24 #include "AliRICHDigit.h"
25 #include "AliRICHRawCluster.h"
26 #include "AliRICHSegmentationV0.h"
27 #include "AliRun.h"
28 #include "TParticle.h"
29 #include "TTree.h"
30 #include "TMath.h"
31 #include "TRandom.h"
32 #include "TH3.h"
33 #include "TH2.h"
34 #include "TCanvas.h"
35 #include <TStyle.h>
36
37
38
39 ClassImp(AliRICHDetect)
40 //___________________________________________
41 AliRICHDetect::AliRICHDetect() : TObject()
42 {
43
44 // Default constructor 
45
46   fc1 = 0;
47   fc2 = 0;
48   fc3 = 0;
49
50 }
51
52 //___________________________________________
53 AliRICHDetect::AliRICHDetect(const char *name, const char *title)
54     : TObject()
55 {
56
57   TStyle *mystyle=new TStyle("Plain","mystyle");
58   mystyle->SetPalette(1,0);
59   mystyle->cd();
60
61
62   fc1= new TCanvas("c1","Reconstructed points",50,50,300,350);
63   fc1->Divide(2,2);
64   fc2= new TCanvas("c2","Reconstructed points after SPOT",370,50,300,350);
65   fc2->Divide(2,2); 
66   fc3= new TCanvas("c3","Used Digits",690,50,300,350);
67   fc4= new TCanvas("c4","Mesh activation data",50,430,600,350);
68   fc4->Divide(2,1);
69  
70
71 }
72
73 //___________________________________________
74 AliRICHDetect::~AliRICHDetect()
75 {
76     
77 // Destructor
78
79 }
80
81
82 void AliRICHDetect::Detect(Int_t nev, Int_t type)
83 {       
84     
85 //
86 // Detection algorithm
87
88
89   //printf("Detection started!\n");
90   Float_t omega,omega1,theta1,phi_relative,steptheta,stepphi,x,y,q=0,z,cx,cy,l,aux1,aux2,aux3,max,radius=0,meanradius=0;
91   Int_t maxi,maxj,maxk;
92   Float_t originalOmega, originalPhi, originalTheta;
93   Float_t binomega, bintheta, binphi;
94   Int_t intomega, inttheta, intphi;
95   Int_t i,j,k;
96
97   AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
98   AliRICHSegmentationV0*  segmentation;
99   AliRICHChamber*       iChamber;
100   AliRICHGeometry*  geometry;
101   
102   iChamber = &(pRICH->Chamber(0));
103   segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
104   geometry=iChamber->GetGeometryModel();
105  
106   
107   //const Float_t Noise_Level=0;                       //Noise Level in percentage of mesh points
108   //const Float_t t=0.6;                               //Softening of Noise Correction (factor)
109   
110   const Float_t kPi=TMath::Pi();                
111   
112   const Float_t kHeight=geometry->GetRadiatorToPads(); //Distance from Radiator to Pads in centimeters
113   //printf("Distance to Pads:%f\n",kHeight);
114  
115   const Int_t kSpot=0;                                 //number of passes with spot algorithm
116   
117   const Int_t kDimensionTheta=100;                     //Matrix dimension for angle Detection
118   const Int_t kDimensionPhi=100;
119   const Int_t kDimensionOmega=100;
120   
121   const Float_t SPOTp=.25;                            //Percentage of spot action
122   const Float_t kMinOmega=.4;
123   const Float_t kMaxOmega=.7;                 //Maximum Cherenkov angle to identify
124   const Float_t kMinTheta=0;
125   const Float_t kMaxTheta=10*kPi/180;   
126   const Float_t kMinPhi=0;
127   const Float_t kMaxPhi=360*kPi/180;
128
129   Float_t rechit[6];                                 //Reconstructed point data
130
131   Int_t ***point = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
132   Int_t ***point1 = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
133   
134   steptheta=(kMaxTheta-kMinTheta)/kDimensionTheta;
135   stepphi=(kMaxPhi-kMinPhi)/kDimensionPhi;
136
137   static TH3F *Points = new TH3F("Points","Reconstructed points 3D",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
138   static TH2F *ThetaPhi = new TH2F("ThetaPhi","Theta-Phi projection",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi);
139   static TH2F *OmegaTheta = new TH2F("OmegaTheta","Omega-Theta projection",kDimensionTheta,0,kDimensionTheta,kDimensionOmega,0,kDimensionOmega);
140   static TH2F *OmegaPhi = new TH2F("OmegaPhi","Omega-Phi projection",kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
141   static TH3F *SpotPoints = new TH3F("Points","Reconstructed points 3D, spot",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
142   static TH2F *SpotThetaPhi = new TH2F("ThetaPhi","Theta-Phi projection, spot",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi);
143   static TH2F *SpotOmegaTheta = new TH2F("OmegaTheta","Omega-Theta projection, spot",kDimensionTheta,0,kDimensionTheta,kDimensionOmega,0,kDimensionOmega);
144   static TH2F *SpotOmegaPhi = new TH2F("OmegaPhi","Omega-Phi projection, spot",kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
145   static TH2F *DigitsXY = new TH2F("DigitsXY","Pads used for reconstruction",150,-25,25,150,-25,25);
146   static TH1F *AngleAct = new TH1F("AngleAct","Activation per angle",100,.45,1);
147   static TH1F *Activation = new TH1F("Activation","Activation per ring",100,0,25);
148   Points->SetXTitle("theta");
149   Points->SetYTitle("phi");
150   Points->SetZTitle("omega");
151   ThetaPhi->SetXTitle("theta");
152   ThetaPhi->SetYTitle("phi");
153   OmegaTheta->SetXTitle("theta");
154   OmegaTheta->SetYTitle("omega");
155   OmegaPhi->SetXTitle("phi");
156   OmegaPhi->SetYTitle("omega");
157   SpotPoints->SetXTitle("theta");
158   SpotPoints->SetYTitle("phi");
159   SpotPoints->SetZTitle("omega");
160   SpotThetaPhi->SetXTitle("theta");
161   SpotThetaPhi->SetYTitle("phi");
162   SpotOmegaTheta->SetXTitle("theta");
163   SpotOmegaTheta->SetYTitle("omega");
164   SpotOmegaPhi->SetXTitle("phi");
165   SpotOmegaPhi->SetYTitle("omega");
166   AngleAct->SetFillColor(5);
167   AngleAct->SetXTitle("rad");
168   AngleAct->SetYTitle("activation");
169   Activation->SetFillColor(5);
170   Activation->SetXTitle("activation");
171
172   Int_t ntracks = (Int_t)pRICH->TreeH()->GetEntries();
173
174   Float_t trackglob[3];
175   Float_t trackloc[3];
176
177   //printf("Area de uma elipse com teta 0 e Omega 45:%f",Area(0,45));
178    
179   Int_t track;
180         
181   for (track=0; track<ntracks;track++) {
182     gAlice->ResetHits();
183     pRICH->TreeH()->GetEvent(track);
184     TClonesArray *pHits  = pRICH->Hits();
185     if (pHits == 0) return;
186     Int_t nhits = pHits->GetEntriesFast();
187     if (nhits == 0) continue;
188     //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
189     AliRICHhit *mHit = 0;
190     //Int_t npoints=0;
191     
192     Int_t counter=0, counter1=0;
193     //Initialization
194     for(i=0;i<kDimensionTheta;i++)
195       {
196         for(j=0;j<kDimensionPhi;j++)
197           {
198             for(k=0;k<kDimensionOmega;k++)
199               {
200                 counter++;
201                 point[i][j][k]=0;
202                 //printf("Dimensions theta:%d, phi:%d, omega:%d",kDimensionTheta,kDimensionPhi,kDimensionOmega);
203                 //printf("Resetting %d %d %d, time %d\n",i,j,k,counter);
204                 //-Noise_Level*(Area(i*kPi/(18*dimension),k*kMaxOmega/dimension)-Area((i-1)*kPi/(18*dimension),(k-1)*kMaxOmega/dimension));
205                 //printf("n-%f",-Noise_Level*(Area(i*kPi/(18*dimension),k*kMaxOmega/dimension)-Area((i-1)*kPi/(18*dimension),(k-1)*kMaxOmega/dimension)));
206               }
207           }
208       }
209
210     Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
211     
212     
213     originalOmega = 0;
214     counter = 0;
215     
216     for (Int_t hit=0;hit<ncerenkovs;hit++) {
217       AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
218       Float_t loss = cHit->fLoss;                 //did it hit the CsI?
219       Float_t production = cHit->fProduction;     //was it produced in freon?  
220       Float_t cherenkov = cHit->fCerenkovAngle;   //production cerenkov angle
221       if (loss == 4 && production == 1)
222         {
223           counter +=1;
224           originalOmega += cherenkov;
225           //printf("%f\n",cherenkov);
226         }
227     }
228     
229     printf("Cerenkovs       : %d\n",counter);
230     
231     if(counter != 0)    //if there are cerenkovs
232       {
233         originalOmega = originalOmega/counter;
234         printf("Original omega: %f\n",originalOmega);
235      
236
237
238         mHit = (AliRICHhit*) pHits->UncheckedAt(0);
239         Int_t nch  = mHit->Chamber();
240         originalTheta = mHit->Theta();
241         originalPhi = mHit->Phi();
242         trackglob[0] = mHit->X();
243         trackglob[1] = mHit->Y();
244         trackglob[2] = mHit->Z();
245         
246         printf("\n--------------------------------------\n");
247         printf("Chamber %d, track %d\n", nch, track);
248
249         
250         iChamber = &(pRICH->Chamber(nch-1));
251     
252         //printf("Nch:%d\n",nch);
253
254         iChamber->GlobaltoLocal(trackglob,trackloc);
255     
256         iChamber->LocaltoGlobal(trackloc,trackglob);
257        
258        
259         cx=trackloc[0];
260         cy=trackloc[2];
261      
262         AliRICHDigit *points = 0;
263         TClonesArray *pDigits = pRICH->DigitsAddress(nch-1);   
264         
265         AliRICHRawCluster *cluster =0;
266         TClonesArray *pClusters = pRICH->RawClustAddress(nch-1);
267
268         Int_t maxcycle=0;
269
270         //digitize from digits
271         
272         if(type==0)
273           {
274             gAlice->TreeD()->GetEvent(0);
275             Int_t ndigits = pDigits->GetEntriesFast();
276             maxcycle=ndigits;
277             printf("Got %d digits\n",ndigits);
278           }
279         
280         //digitize from clusters
281         
282         if(type==1)
283           {
284             Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
285             gAlice->TreeR()->GetEvent(nent-1);
286             Int_t nclusters = pClusters->GetEntriesFast();
287             maxcycle=nclusters;
288             printf("Got %d clusters\n",nclusters);
289           }
290
291
292
293         
294         counter=0;
295         printf("Starting calculations\n");
296         printf("           Start                                                                                              Finish\n");
297         printf("Progress:     ");
298         for(Float_t theta=0;theta<kMaxTheta;theta+=steptheta)
299           {             
300             printf(".");
301             for(Float_t phi=0;phi<=kMaxPhi;phi+=stepphi)
302               {         
303                 //printf("Phi:%3.1f\n", phi*180/kPi);
304                 counter1=0;
305                 for (Int_t cycle=0;cycle<maxcycle;cycle++)
306                   {     
307                     
308                     if(type==0)
309                       {
310                         points=(AliRICHDigit*) pDigits->UncheckedAt(cycle);
311                         segmentation->GetPadC(points->PadX(), points->PadY(),x, y, z);
312                       }
313                     
314                     if(type==1)
315                       {
316                         cluster=(AliRICHRawCluster*) pClusters->UncheckedAt(cycle);
317                         x=cluster->fX;
318                         y=cluster->fY;
319                         q=cluster->fQ;
320                       }
321                     
322                     if(type ==0 || q > 100)
323                       
324                       {
325                         
326                         x=x-cx;
327                         y=y-cy;
328                         radius=TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
329                         
330                         //calculation of relative phi angle of digit
331                         
332                         phi_relative = acos(y/radius);
333                         phi_relative = TMath::Abs(phi_relative - phi);
334                         
335                         
336                         if(radius>4)
337                           {
338                             meanradius+=radius;
339                             counter++;
340                             
341                             
342                             //if (radius < (2*kHeight*tan(theta+kMaxOmega)))
343                             if (radius < (2*kHeight*tan(kMaxOmega)))
344                               //if(Fiducial(x,y,theta,phi,kHeight,kMaxOmega,kMinOmega))
345                               {
346                                 
347                                 if(phi==0)
348                                   {
349                                 //printf("Radius: %f, Max Radius: %f\n",radius,2*kHeight*tan(theta+kMaxOmega)*3/4);
350                                 //printf("Loaded digit %d with coordinates x:%f, y%f\n",dig,x,y);
351                                 //printf("Using digit %d, for theta %f\n",dig,theta);
352                                   }
353                                 
354                                 counter1++;
355                                 
356                                 l=kHeight/cos(theta);
357                                 
358                                 //main calculation
359                                 
360                                 DigitsXY->Fill(x,y,(float) 1);
361                                 
362                                 theta1=SnellAngle(theta);
363                                 
364                                 aux1=-y*sin(phi)+x*cos(phi);
365                                 aux2=y*cos(phi)+x*sin(phi);
366                                 aux3=( TMath::Power(aux1,2)+TMath::Power(cos(theta1)*aux2 ,2))/TMath::Power(sin(theta1)*aux2+l,2);
367                                 omega=atan(sqrt(aux3));
368                                 
369                                 //omega is distorted, theta1 is distorted
370                                 
371                                 if(InvSnellAngle(omega+TMath::Abs(cos(phi_relative))*theta1)<999)
372                                   {
373                                     omega1=InvSnellAngle(omega+TMath::Abs(cos(phi_relative))*theta);
374                                     theta1=InvSnellAngle(theta1);
375                                     
376                                   }
377                                 else
378                                   {
379                                     omega1=0;
380                                     theta1=0;
381                                   }
382                                 
383                                 if(omega1<kMaxOmega && omega1>kMinOmega)
384                                   {
385                                 //printf("Omega found:%f\n",omega);
386                                     omega1=omega1-kMinOmega;
387                                     
388                                 //printf("Omega: %f Theta: %3.1f Phi:%3.1f\n",omega, theta*180/kPi, phi*180/kPi);
389                                     
390                                     bintheta=theta1*kDimensionTheta/kMaxTheta;
391                                     binphi=phi*kDimensionPhi/kMaxPhi;
392                                     binomega=omega1*kDimensionOmega/(kMaxOmega-kMinOmega);
393                                     
394                                     if(Int_t(bintheta+0.5)==Int_t(bintheta))
395                                       inttheta=Int_t(bintheta);
396                                     else
397                                       inttheta=Int_t(bintheta+0.5);
398                                     
399                                     if(Int_t(binomega+0.5)==Int_t(binomega))
400                                       intomega=Int_t(binomega);
401                                     else
402                                       intomega=Int_t(binomega+0.5);
403                                     
404                                     if(Int_t(binphi+0.5)==Int_t(binphi))
405                                       intphi=Int_t(binphi);
406                                     else
407                                       intphi=Int_t(binphi+0.5);
408                                     
409                                 //printf("Point added at %d %d %d\n",inttheta,intphi,intomega);
410                                     
411                                     if(type==0)
412                                       point[inttheta][intphi][intomega]+=1;
413                                     if(type==1)
414                                       point[inttheta][intphi][intomega]+=(Int_t)(q);
415                                     
416                                 //printf("Omega stored:%d\n",intomega);
417                                     Points->Fill(inttheta,intphi,intomega,(float) 1);
418                                     ThetaPhi->Fill(inttheta,intphi,(float) 1);
419                                     OmegaTheta->Fill(inttheta,intomega,(float) 1);
420                                     OmegaPhi->Fill(intphi,intomega,(float) 1);
421                                 //printf("Filling at %d %d %d\n",Int_t(theta*kDimensionTheta/kMaxTheta),Int_t(phi*kDimensionPhi/kMaxPhi),Int_t(omega*kDimensionOmega/kMaxOmega));
422                                   }
423                                 //if(omega<kMaxOmega)point[Int_t(theta)][Int_t(phi)][Int_t(omega)]+=1;
424                               }
425                           }
426                       }
427                     
428                   }
429                 //printf("Used %d digits for theta %3.1f\n",counter1, theta*180/kPi);
430               }
431             
432           }
433
434         printf("\n");
435
436         meanradius=meanradius/counter;
437         //printf("Mean radius:%f, counter:%d\n",meanradius,counter);
438         rechit[5]=meanradius;
439         //printf("Used %d digits\n",counter1);
440         //printf("\n");
441
442         if(nev<2)
443           {
444             if(nev==0)
445               {
446                 fc1->cd(1);
447                 Points->Draw("colz");
448                 fc1->cd(2);
449                 ThetaPhi->Draw("colz");
450                 fc1->cd(3);
451                 OmegaTheta->Draw("colz");
452                 fc1->cd(4);
453                 OmegaPhi->Draw("colz");
454                 fc3->cd();
455                 DigitsXY->Draw("colz");
456               }
457             else
458               {
459                 //fc1->cd(1);
460                 //Points->Draw("same");
461                 //fc1->cd(2);
462                 //ThetaPhi->Draw("same");
463                 //fc1->cd(3);
464                 //OmegaTheta->Draw("same");
465                 //fc1->cd(4);
466                 //OmegaPhi->Draw("same");       
467               }
468           }
469         
470     
471         //SPOT execute twice
472         for(Int_t s=0;s<kSpot;s++)
473           {
474             printf("     Applying Spot algorithm, pass %d\n", s);
475         
476         //buffer copy
477         for(i=0;i<=kDimensionTheta;i++)
478           {
479             for(j=0;j<=kDimensionPhi;j++)
480               {
481                 for(k=0;k<=kDimensionOmega;k++)
482                   {
483                     point1[i][j][k]=point[i][j][k];     
484                   }
485               }
486           }
487
488         //SPOT algorithm                        
489         for(i=1;i<kDimensionTheta-1;i++)
490           {
491             for(j=1;j<kDimensionPhi-1;j++)
492               {
493                 for(k=1;k<kDimensionOmega-1;k++)
494                   {
495                     if((point[i][k][j]>point[i-1][k][j])&&(point[i][k][j]>point[i+1][k][j])&&
496                        (point[i][k][j]>point[i][k-1][j])&&(point[i][k][j]>point[i][k+1][j])&&
497                        (point[i][k][j]>point[i][k][j-1])&&(point[i][k][j]>point[i][k][j+1]))
498                       {
499                         //cout<<"SPOT"<<endl;
500                         //Execute SPOT on point                                                                                         
501                         point1[i][j][k]+=Int_t(SPOTp*(point[i-1][k][j]+point[i+1][k][j]+point[i][k-1][j]+point[i][k+1][j]+point[i][k][j-1]+point[i][k][j+1]));    
502                         point1[i-1][k][j]=Int_t(SPOTp*point[i-1][k][j]);
503                         point1[i+1][k][j]=Int_t(SPOTp*point[i+1][k][j]);
504                         point1[i][k-1][j]=Int_t(SPOTp*point[i][k-1][j]);
505                         point1[i][k+1][j]=Int_t(SPOTp*point[i][k+1][j]);
506                         point1[i][k][j-1]=Int_t(SPOTp*point[i][k][j-1]);
507                         point1[i][k][j+1]=Int_t(SPOTp*point[i][k][j+1]);
508                       }
509                   }
510               }
511           }
512         
513         //copy from buffer copy
514         counter1=0;
515         for(i=1;i<kDimensionTheta;i++)
516           {
517             for(j=1;j<kDimensionPhi;j++)
518               {
519                 for(k=1;k<kDimensionOmega;k++)
520                   {
521                     point[i][j][k]=point1[i][j][k];
522                     if(nev<20)
523                       {
524                         if(s==kSpot-1)
525                           {
526                             if(point1[i][j][k] != 0)
527                               {
528                                 SpotPoints->Fill(i,j,k,(float) point1[i][j][k]);
529                                 //printf("Random number %f\n",random->Rndm(2));
530                                 //if(random->Rndm() < .2)
531                                   //{
532                                 SpotThetaPhi->Fill(i,j,(float) point1[i][j][k]);
533                                 SpotOmegaTheta->Fill(i,k,(float) point1[i][j][k]);
534                                 SpotOmegaPhi->Fill(j,k,(float) point1[i][j][k]);
535                                     counter1++;
536                                   //}
537                                 //printf("Filling at %d %d %d value %f\n",i,j,k,(float) point1[i][j][k]);
538                               }
539                           }
540                       }
541                     //if(point1[i][j][k] != 0)
542                       //printf("Last transfer point: %d, point1, %d\n",point[i][j][k],point1[i][j][k]);
543                   }
544               }
545           }
546       }
547     
548     //printf("Filled %d cells\n",counter1);
549
550         if(nev<2)
551           {
552             if(nev==0)
553               {
554                 fc2->cd(1);
555                 SpotPoints->Draw("colz");
556                 fc2->cd(2);
557                 SpotThetaPhi->Draw("colz");
558                 fc2->cd(3);
559                 SpotOmegaTheta->Draw("colz");
560                 fc2->cd(4);
561                 SpotOmegaPhi->Draw("colz");
562               }
563             else
564               {
565                 //fc2->cd(1);
566                 //SpotPoints->Draw("same");
567                 //fc2->cd(2);
568                 //SpotThetaPhi->Draw("same");
569                 //fc2->cd(3);
570                 //SpotOmegaTheta->Draw("same");
571                 //fc2->cd(4);
572                 //SpotOmegaPhi->Draw("same");   
573               }
574           }
575     
576     
577         //Identification is equivalent to maximum determination
578         max=0;maxi=0;maxj=0;maxk=0;
579         
580         printf("     Proceeding to identification");
581         
582         for(i=0;i<kDimensionTheta;i++)
583           for(j=0;j<kDimensionPhi;j++)
584             for(k=0;k<kDimensionOmega;k++)
585               if(point[i][j][k]>max)
586                 {
587                   //cout<<"maxi="<<i*90/dimension<<" maxj="<<j*90/dimension<<" maxk="<<k*kMaxOmega/dimension*180/kPi<<" max="<<max<<endl;
588                 maxi=i;maxj=j;maxk=k;
589                 max=point[i][j][k];
590                 printf(".");
591                 //printf("Max Omega %d, Max Theta %d, Max Phi %d (%d counts)\n",maxk,maxi,maxj,max);
592                 }
593         printf("\n");
594     
595         Float_t FinalOmega = maxk*(kMaxOmega-kMinOmega)/kDimensionOmega; 
596         Float_t FinalTheta = maxi*kMaxTheta/kDimensionTheta;
597         Float_t FinalPhi = maxj*kMaxPhi/kDimensionPhi;
598
599         FinalOmega += kMinOmega;
600     
601         //printf("Detected angle for height %3.1f and for center %3.1f %3.1f:%f\n",h,cx,cy,maxk*kPi/(kDimensionTheta*4));
602         printf("     Indentified angles: cerenkov - %f, theta - %3.1f, phi - %3.1f (%f activation)\n", FinalOmega, FinalTheta*180/kPi, FinalPhi*180/kPi, max);
603         //printf("Detected angle for height %3.1f and for center %3.1f %3.1f:%f\n",kHeight,cx,cy,maxk);
604
605         AngleAct->Fill(FinalOmega, (float) max);
606         Activation->Fill(max, (float) 1);
607
608         if(nev==0)
609               {
610                 fc4->cd(1);
611                 AngleAct->Draw();
612                 fc4->cd(2);
613                 Activation->Draw();
614               }
615             else
616               {
617                 fc4->cd(1);
618                 AngleAct->Draw("same");
619                 fc4->cd(2);
620                 Activation->Draw("same");
621               }
622           
623
624         //fscanf(omegas,"%f",&realomega);
625         //fscanf(thetas,"%f",&realtheta);
626         //printf("Real Omega: %f",realomega);                   
627         //cout<<"Detected:theta="<<maxi*90/kDimensionTheta<<"phi="<<maxj*90/kDimensionPhi<<"omega="<<maxk*kMaxOmega/kDimensionOmega*180/kPi<<" OmegaError="<<fabs(maxk*kMaxOmega/kDimensionOmega*180/kPi-realomega)<<" ThetaError="<<fabs(maxi*90/kDimensionTheta-realtheta)<<endl<<endl;               
628     
629         //fprintf(results,"Center Coordinates, cx=%6.2f cy=%6.2f, Real Omega=%6.2f, Detected Omega=%6.2f, Omega Error=%6.2f Theta Error=%6.2f\n",cx,cy,realomega,maxk*kMaxOmega/kDimensionOmega*180/kPi,fabs(maxk*kMaxOmega/kDimensionOmega*180/kPi-realomega),fabs(maxi*90/kDimensionTheta-realtheta));
630     
631     /*for(j=0;j<np;j++)
632       pointpp(maxj*90/kDimensionTheta,maxi*90/kDimensionPhi,maxk*kMaxOmega/kDimensionOmega*180/kPi,cx,cy);//Generates a point on the elipse*/               
633
634
635     //Start filling rec. hits
636     
637     rechit[0] = FinalTheta;
638     rechit[1] = FinalPhi - 90*kPi/180;
639     rechit[2] = FinalOmega;
640     rechit[3] = cx;
641     rechit[4] = cy;
642
643     //CreatePoints(FinalTheta, 270*kPi/180 + FinalPhi, FinalOmega, kHeight);
644        
645     printf ("track %d, theta %f, phi %f, omega %f\n\n\n",track,rechit[0],rechit[1]*180/kPi,rechit[2]);
646     
647     // fill rechits
648     pRICH->AddRecHit3D(nch-1,rechit,originalOmega, originalTheta, originalPhi);
649     //printf("rechit %f %f %f %f %f\n",rechit[0],rechit[1],rechit[2],rechit[3],rechit[4]);
650     //printf("Chamber:%d",nch);
651       }
652
653     else   //if no cerenkovs
654       
655       {
656         
657         rechit[0] = 0;
658         rechit[1] = 0;
659         rechit[2] = 0;
660         rechit[3] = 0;
661         rechit[4] = 0;
662
663       }
664
665   }
666
667   if(type==1)  //reco from clusters
668     {
669       pRICH->ResetRawClusters();
670       //Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
671       //gAlice->TreeR()->GetEvent(track);
672       //printf("Going to branch %d\n",track);
673       //gAlice->GetEvent(nev);
674     }
675
676         
677     //printf("\n\n\n\n");
678     gAlice->TreeR()->Fill();
679   TClonesArray *fRec;
680   for (i=0;i<kNCH;i++) {
681     fRec=pRICH->RecHitsAddress3D(i);
682     int ndig=fRec->GetEntriesFast();
683     printf ("Chamber %d, rings %d\n",i+1,ndig);
684   }
685   pRICH->ResetRecHits3D();
686
687   free_i3tensor(point,0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
688   free_i3tensor(point1,0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
689 }
690
691
692
693 Float_t AliRICHDetect:: Area(Float_t theta,Float_t omega)
694 {
695
696 //
697 // Calculates area of an ellipse for given incidence angles    
698
699
700     Float_t area;
701     const Float_t kHeight=9.25;                       //Distance from Radiator to Pads in pads
702     
703     area=TMath::Pi()*TMath::Power(kHeight*tan(omega),2)/TMath::Power(TMath::Power(cos(theta),2)-TMath::Power(tan(omega)*sin(theta),2),3/2);
704     
705     return (area);
706 }
707
708
709 Int_t ***AliRICHDetect::i3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
710 // allocate a Int_t 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] 
711 {
712   long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
713   Int_t ***t;
714   
715   int NR_END=1; 
716
717   // allocate pointers to pointers to rows 
718   t=(Int_t ***) malloc((size_t)((nrow+NR_END)*sizeof(Int_t**)));
719   if (!t) printf("allocation failure 1 in f3tensor()");
720   t += NR_END;
721   t -= nrl;
722   
723   // allocate pointers to rows and set pointers to them 
724   t[nrl]=(Int_t **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(Int_t*)));
725   if (!t[nrl]) printf("allocation failure 2 in f3tensor()");
726   t[nrl] += NR_END;
727   t[nrl] -= ncl;
728   
729   // allocate rows and set pointers to them 
730   t[nrl][ncl]=(Int_t *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(Int_t)));
731   if (!t[nrl][ncl]) printf("allocation failure 3 in f3tensor()");
732   t[nrl][ncl] += NR_END;
733   t[nrl][ncl] -= ndl;
734   
735   for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
736   for(i=nrl+1;i<=nrh;i++) {
737     t[i]=t[i-1]+ncol;
738     t[i][ncl]=t[i-1][ncl]+ncol*ndep;
739     for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
740   }
741   
742   // return pointer to array of pointers to rows 
743   return t;
744 }
745
746 void AliRICHDetect::free_i3tensor(int ***t, long nrl, long nrh, long ncl, long nch,long ndl, long ndh)
747 // free a Int_t f3tensor allocated by i3tensor()
748 {
749   int NR_END=1; 
750
751   free((char*) (t[nrl][ncl]+ndl-NR_END));
752   free((char*) (t[nrl]+ncl-NR_END));
753   free((char*) (t+nrl-NR_END));
754 }
755
756
757 Float_t AliRICHDetect:: SnellAngle(Float_t iangle)
758
759
760 // Compute the Snell angle
761
762   Float_t nfreon  = 1.295;
763   Float_t nquartz = 1.585;
764   Float_t ngas = 1;
765
766   Float_t sinrangle;
767   Float_t rangle;
768   Float_t a1, a2;
769
770   a1=nfreon/nquartz;
771   a2=nquartz/ngas;
772
773   sinrangle = a1*a2*sin(iangle);
774
775   if(sinrangle>1.) {
776      rangle = 999.;
777      return rangle;
778   }
779   
780   rangle = asin(sinrangle);  
781   return rangle;
782 }
783
784 Float_t AliRICHDetect:: InvSnellAngle(Float_t rangle)
785
786
787 // Compute the inverse Snell angle
788
789   Float_t nfreon  = 1.295;
790   Float_t nquartz = 1.585;
791   Float_t ngas = 1;
792
793   Float_t siniangle;
794   Float_t iangle;
795   Float_t a1,a2;
796
797   a1=nfreon/nquartz;
798   a2=nquartz/ngas;
799
800   siniangle = sin(rangle)/(a1*a2);
801   iangle = asin(siniangle);
802
803   if(siniangle>1.) {
804      iangle = 999.;
805      return iangle;
806   }
807   
808   iangle = asin(siniangle);
809   return iangle;
810 }
811
812
813
814 //________________________________________________________________________________
815 void AliRICHDetect::CreatePoints(Float_t theta, Float_t phi, Float_t omega, Float_t h)
816 {
817   
818   // Create points along the ellipse equation
819   
820   Int_t s1,s2;
821   Float_t fiducial=h*TMath::Tan(omega+theta), l=h/TMath::Cos(theta), xtrial, y=0, c0, c1, c2;
822   //TRandom *random=new TRandom();
823   
824   static TH2F *REllipse = new TH2F("REllipse","Reconstructed ellipses",150,-25,25,150,-25,25);
825
826   for(Float_t i=0;i<1000;i++)
827     {
828       
829       Float_t counter=0;
830       
831       c0=0;c1=0;c2=0;
832       while((c1*c1-4*c2*c0)<=0 && counter<1000)
833         {
834           //Choose which side to go...
835           if(i>250 && i<750) s1=1; 
836           //if (gRandom->Rndm(1)>.5) s1=1;
837           else s1=-1;
838           //printf("s1:%d\n",s1);
839           //Trial a y
840           y=s1*i*gRandom->Rndm(Int_t(fiducial/50));
841           //printf("Fiducial %f  for omega:%f theta:%f phi:%f\n",fiducial,omega,theta,fphi);
842           Float_t alfa1=theta;
843           Float_t theta1=phi;
844           Float_t omega1=omega;
845           
846           //Solve the eq for a trial x
847           c0=-TMath::Power(y*TMath::Cos(alfa1)*TMath::Cos(theta1),2)-TMath::Power(y*TMath::Sin(alfa1),2)+TMath::Power(l*TMath::Tan(omega1),2)+2*l*y*TMath::Cos(alfa1)*TMath::Sin(theta1)*TMath::Power(TMath::Tan(omega1),2)+TMath::Power(y*TMath::Cos(alfa1)*TMath::Sin(theta1)*TMath::Tan(omega1),2);
848           c1=2*y*TMath::Cos(alfa1)*TMath::Sin(alfa1)-2*y*TMath::Cos(alfa1)*TMath::Power(TMath::Cos(theta1),2)*TMath::Sin(alfa1)+2*l*TMath::Sin(alfa1)*TMath::Sin(theta1)*TMath::Power(TMath::Tan(omega1),2)+2*y*TMath::Cos(alfa1)*TMath::Sin(alfa1)*TMath::Power(TMath::Sin(theta1),2)*TMath::Power(TMath::Tan(omega1),2);
849           c2=-TMath::Power(TMath::Cos(alfa1),2)-TMath::Power(TMath::Cos(theta1)*TMath::Sin(alfa1),2)+TMath::Power(TMath::Sin(alfa1)*TMath::Sin(theta1)*TMath::Tan(omega1),2);
850           //cout<<"Trial: y="<<y<<"c0="<<c0<<" c1="<<c1<<" c2="<<c2<<endl;
851           //printf("Result:%f\n\n",c1*c1-4*c2*c0);
852           //i+=.01;
853           counter +=1;
854         }
855       
856       if (counter>=1000)
857         y=0; 
858
859       //Choose which side to go...
860       //if(gRandom->Rndm(1)>.5) s=1; 
861       //else s=-1;
862       if(i>500) s2=1;
863       //if (gRandom->Rndm(1)>.5) s2=1;
864       else s2=-1;
865       xtrial=(-c1+s2*TMath::Sqrt(c1*c1-4*c2*c0))/(2*c2);
866       //cout<<"x="<<xtrial<<" y="<<cy+y<<endl;
867       //printf("Coordinates: %f %f\n",xtrial,fCy+y);
868
869       REllipse->Fill(xtrial,y);
870
871       //printf("Coordinates: %f %f %f\n",vectorGlob[0],vectorGlob[1],vectorGlob[2]);
872     }
873   
874   fc3->cd(2);
875   REllipse->Draw();
876 }
877
878 Int_t AliRICHDetect::Fiducial(Float_t rotx, Float_t roty, Float_t theta, Float_t phi, Float_t height, Float_t maxOmega, Float_t minOmega)
879 {
880
881   Float_t x,y,y1,h,omega1,omega2;
882
883   Float_t a,b, offset;
884   Float_t a1,b1, offset1;
885
886   h = height;
887
888   //refraction calculation
889   
890   theta = SnellAngle(theta);
891   phi = phi - TMath::Pi();
892   omega1 = SnellAngle(maxOmega);
893   omega2 = SnellAngle(minOmega);
894
895   //maximum zone
896   a = h*(tan(omega1+theta)+tan(omega1-theta))/2;
897   b = h*tan(omega1);
898  
899   offset = h*(tan(omega1+theta)-tan(omega1-theta))/2;
900   
901  
902   //minimum zone
903   a1 = h*(tan(omega2+theta)+tan(omega2-theta))/2;
904   b1 = h*tan(omega2);
905  
906   offset1 = h*(tan(omega2+theta)-tan(omega2-theta))/2;
907   
908
909   //rotation to phi=0
910   x = rotx*cos(phi)+roty*sin(phi);
911   y = -rotx*sin(phi)+roty*cos(phi) - offset;
912   y1 = -rotx*sin(phi)+roty*cos(phi) - offset1;
913
914   
915   if(x*x/a+y*y/b<1 && x*x/a1+y1*y1/b1>1)
916     return 1;
917   else
918     return 0;
919
920 }