Second reconstruction algorithm. (J. Barbosa)
[u/mrichter/AliRoot.git] / RICH / AliRICHDetectV1.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 /*
17   $Log$
18  
19 */
20
21 #include <stdlib.h>
22
23 #include "AliRICHDetectV1.h"
24 #include "AliRICH.h"
25 #include "AliRICHPoints.h"
26 #include "AliRICHDetect.h"
27 #include "AliRICHHit.h"
28 #include "AliRICHDigit.h"
29 #include "AliRICHRawCluster.h"
30 #include "AliRICHCerenkov.h"
31 #include "AliRICHSegmentationV0.h"
32 #include "AliRun.h"
33 #include "TParticle.h"
34 #include "TTree.h"
35 #include "TMath.h"
36 #include "TRandom.h"
37 #include "TH3.h"
38 #include "TH2.h"
39 #include "TCanvas.h"
40 #include <TStyle.h>
41
42
43 #include "malloc.h"
44
45
46 ClassImp(AliRICHDetectV1)
47
48
49 //___________________________________________
50 AliRICHDetectV1::AliRICHDetectV1() : AliRICHDetect()
51 {
52
53 // Default constructor 
54
55   fc1 = 0;
56   fc2 = 0;
57   fc3 = 0;
58
59 }
60
61 //___________________________________________
62 AliRICHDetectV1::AliRICHDetectV1(const char *name, const char *title)
63     : AliRICHDetect()
64 {
65
66   TStyle *mystyle=new TStyle("Plain","mystyle");
67   mystyle->SetPalette(1,0);
68   mystyle->cd();
69
70
71   fc1= new TCanvas("c1","Reconstructed points",50,50,300,350);
72   fc1->Divide(2,2);
73   fc2= new TCanvas("c2","Reconstructed points after SPOT",370,50,300,350);
74   fc2->Divide(2,2); 
75   fc3= new TCanvas("c3","Used Digits",690,50,300,350);
76   fc4= new TCanvas("c4","Mesh activation data",50,430,600,350);
77   fc4->Divide(2,1);
78  
79
80 }
81
82 //___________________________________________
83 AliRICHDetectV1::~AliRICHDetectV1()
84 {
85     
86 // Destructor
87
88 }
89
90
91 void AliRICHDetectV1::Detect(Int_t nev, Int_t type)
92 {       
93     
94 //
95 // Detection algorithm
96
97
98   //printf("Detection started!\n");
99   Float_t omega,omega1,theta1,x,y,q=0,z,cx,cy,max,radius=0,meanradius=0;
100   Int_t maxi,maxj,maxk;
101   Float_t originalOmega, originalPhi, originalTheta;
102   Float_t steptheta,stepphi,stepomega;
103   Float_t binomega, bintheta, binphi;
104   Int_t intomega, inttheta, intphi;
105   Float_t maxRadius,minRadius,eccentricity,angularadius,offset,phi_relative;
106   Int_t i,j,k;
107
108   AliRICH *pRICH  = (AliRICH*)gAlice->GetDetector("RICH");
109   AliRICHSegmentationV0*  segmentation;
110   AliRICHChamber*       iChamber;
111   AliRICHGeometry*  geometry;
112   
113   iChamber = &(pRICH->Chamber(0));
114   segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel(0);
115   geometry=iChamber->GetGeometryModel();
116  
117   
118   //const Float_t Noise_Level=0;                       //Noise Level in percentage of mesh points
119   //const Float_t t=0.6;                               //Softening of Noise Correction (factor)
120   
121   const Float_t kPi=TMath::Pi();                
122   
123   const Float_t kHeight=geometry->GetRadiatorToPads(); //Distance from Radiator to Pads in centimeters
124   //printf("Distance to Pads:%f\n",kHeight);
125  
126   const Int_t kSpot=2;                                 //number of passes with spot algorithm
127   const Int_t activ_tresh=0;                           //activation treshold to identify a track
128   
129   const Int_t kDimensionTheta=2;                       //Matrix dimension for angle Detection
130   const Int_t kDimensionPhi=2;
131   const Int_t kDimensionOmega=50;
132   
133   const Float_t SPOTp=.25;                            //Percentage of spot action
134   const Float_t kMinOmega=.6;
135   const Float_t kMaxOmega=.7;                 //Maximum Cherenkov angle to identify
136   const Float_t kMinTheta=0;
137   const Float_t kMaxTheta=0.5*kPi/180;  
138   const Float_t kMinPhi=0;
139   const Float_t kMaxPhi=20*kPi/180;
140
141   const Float_t sigma=0.5;                          //half thickness of fiducial band in cm
142
143   Float_t rechit[6];                                 //Reconstructed point data
144
145   Int_t ***point = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
146   Int_t ***point1 = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
147   
148   steptheta=(kMaxTheta-kMinTheta)/kDimensionTheta;
149   stepphi=(kMaxPhi-kMinPhi)/kDimensionPhi;
150   stepomega=(kMaxOmega-kMinOmega)/kDimensionOmega;
151
152   static TH3F *Points = new TH3F("Points","Reconstructed points 3D",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
153   static TH2F *ThetaPhi = new TH2F("ThetaPhi","Theta-Phi projection",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi);
154   static TH2F *OmegaTheta = new TH2F("OmegaTheta","Omega-Theta projection",kDimensionTheta,0,kDimensionTheta,kDimensionOmega,0,kDimensionOmega);
155   static TH2F *OmegaPhi = new TH2F("OmegaPhi","Omega-Phi projection",kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
156   static TH3F *SpotPoints = new TH3F("Points","Reconstructed points 3D, spot",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
157   static TH2F *SpotThetaPhi = new TH2F("ThetaPhi","Theta-Phi projection, spot",kDimensionTheta,0,kDimensionTheta,kDimensionPhi,0,kDimensionPhi);
158   static TH2F *SpotOmegaTheta = new TH2F("OmegaTheta","Omega-Theta projection, spot",kDimensionTheta,0,kDimensionTheta,kDimensionOmega,0,kDimensionOmega);
159   static TH2F *SpotOmegaPhi = new TH2F("OmegaPhi","Omega-Phi projection, spot",kDimensionPhi,0,kDimensionPhi,kDimensionOmega,0,kDimensionOmega);
160   static TH2F *DigitsXY = new TH2F("DigitsXY","Pads used for reconstruction",150,-25,25,150,-25,25);
161   static TH1F *AngleAct = new TH1F("AngleAct","Activation per angle",100,.45,1);
162   static TH1F *Activation = new TH1F("Activation","Activation per ring",100,0,25);
163   Points->SetXTitle("theta");
164   Points->SetYTitle("phi");
165   Points->SetZTitle("omega");
166   ThetaPhi->SetXTitle("theta");
167   ThetaPhi->SetYTitle("phi");
168   OmegaTheta->SetXTitle("theta");
169   OmegaTheta->SetYTitle("omega");
170   OmegaPhi->SetXTitle("phi");
171   OmegaPhi->SetYTitle("omega");
172   SpotPoints->SetXTitle("theta");
173   SpotPoints->SetYTitle("phi");
174   SpotPoints->SetZTitle("omega");
175   SpotThetaPhi->SetXTitle("theta");
176   SpotThetaPhi->SetYTitle("phi");
177   SpotOmegaTheta->SetXTitle("theta");
178   SpotOmegaTheta->SetYTitle("omega");
179   SpotOmegaPhi->SetXTitle("phi");
180   SpotOmegaPhi->SetYTitle("omega");
181   AngleAct->SetFillColor(5);
182   AngleAct->SetXTitle("rad");
183   AngleAct->SetYTitle("activation");
184   Activation->SetFillColor(5);
185   Activation->SetXTitle("activation");
186
187   Int_t ntracks = (Int_t)gAlice->TreeH()->GetEntries();
188    
189   Float_t trackglob[3];
190   Float_t trackloc[3];
191
192   //printf("Area de uma elipse com teta 0 e Omega 45:%f",Area(0,45));
193    
194   Int_t track;
195         
196   for (track=0; track<ntracks;track++) {
197     gAlice->ResetHits();
198     gAlice->TreeH()->GetEvent(track);
199     TClonesArray *pHits  = pRICH->Hits();
200     if (pHits == 0) return;
201     Int_t nhits = pHits->GetEntriesFast();
202     if (nhits == 0) continue;
203     //Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
204     AliRICHHit *mHit = 0;
205     //Int_t npoints=0;
206     
207     Int_t counter=0, counter1=0;
208     //Initialization
209     for(i=0;i<kDimensionTheta;i++)
210       {
211         for(j=0;j<kDimensionPhi;j++)
212           {
213             for(k=0;k<kDimensionOmega;k++)
214               {
215                 counter++;
216                 point[i][j][k]=0;
217                 //printf("Dimensions theta:%d, phi:%d, omega:%d",kDimensionTheta,kDimensionPhi,kDimensionOmega);
218                 //printf("Resetting %d %d %d, time %d\n",i,j,k,counter);
219                 //-Noise_Level*(Area(i*kPi/(18*dimension),k*kMaxOmega/dimension)-Area((i-1)*kPi/(18*dimension),(k-1)*kMaxOmega/dimension));
220                 //printf("n-%f",-Noise_Level*(Area(i*kPi/(18*dimension),k*kMaxOmega/dimension)-Area((i-1)*kPi/(18*dimension),(k-1)*kMaxOmega/dimension)));
221               }
222           }
223       }
224
225     Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
226     
227     
228     originalOmega = 0;
229     counter = 0;
230     
231     for (Int_t hit=0;hit<ncerenkovs;hit++) {
232       AliRICHCerenkov* cHit = (AliRICHCerenkov*) pRICH->Cerenkovs()->UncheckedAt(hit);
233       Float_t loss = cHit->fLoss;                 //did it hit the CsI?
234       Float_t production = cHit->fProduction;     //was it produced in freon?  
235       Float_t cherenkov = cHit->fCerenkovAngle;   //production cerenkov angle
236       if (loss == 4 && production == 1)
237         {
238           counter +=1;
239           originalOmega += cherenkov;
240           //printf("%f\n",cherenkov);
241         }
242     }
243
244     originalOmega = originalOmega/counter;
245     
246     //printf("Cerenkovs       : %d\n",counter);
247     
248     mHit = (AliRICHHit*) pHits->UncheckedAt(0);
249     Int_t nch  = mHit->Chamber();
250     originalTheta = mHit->Theta();
251     originalPhi = mHit->Phi();
252     trackglob[0] = mHit->X();
253     trackglob[1] = mHit->Y();
254     trackglob[2] = mHit->Z();
255     
256     
257     printf("\n--------------------------------------\n");
258     printf("Chamber %d, track %d\n", nch, track);
259     printf("Original omega: %f\n",originalOmega);
260     
261     iChamber = &(pRICH->Chamber(nch-1));
262     
263     printf("Nch:%d x:%f y:%f\n",nch,trackglob[0],trackglob[2]);
264     
265     iChamber->GlobaltoLocal(trackglob,trackloc);
266     
267     iChamber->LocaltoGlobal(trackloc,trackglob);
268     
269     
270     cx=trackloc[0];
271     cy=trackloc[2];
272
273     printf("cy:%f  ", cy);
274     
275     if(counter != 0)    //if there are cerenkovs
276       {
277         
278         AliRICHDigit *points = 0;
279         TClonesArray *pDigits = pRICH->DigitsAddress(nch-1);   
280         
281         AliRICHRawCluster *cluster =0;
282         TClonesArray *pClusters = pRICH->RawClustAddress(nch-1);
283
284         Int_t maxcycle=0;
285
286         //digitize from digits
287         
288         if(type==0)
289           {
290             gAlice->TreeD()->GetEvent(0);
291             Int_t ndigits = pDigits->GetEntriesFast();
292             maxcycle=ndigits;
293             //printf("Got %d digits\n",ndigits);
294           }
295         
296         //digitize from clusters
297         
298         if(type==1)
299           {
300             Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
301             gAlice->TreeR()->GetEvent(nent-1);
302             Int_t nclusters = pClusters->GetEntriesFast();
303             maxcycle=nclusters;
304             //printf("Got %d clusters\n",nclusters);
305           }
306
307
308
309         
310         counter=0;
311         printf("Starting calculations\n");
312         printf("           Start                                                                                              Finish\n");
313         printf("Progress:     ");
314
315
316         for(Float_t theta=0;theta<kMaxTheta;theta+=steptheta)
317           {             
318             printf(".");
319             for(Float_t phi=0;phi<=kMaxPhi;phi+=stepphi)
320               {         
321                 for(omega=kMinOmega;omega<=kMaxOmega;omega+=stepomega)
322                   {             
323                     //printf("Entering angle cycle\n");
324                     omega1=SnellAngle(omega);
325                     theta1=SnellAngle(theta);
326                     
327                     maxRadius = kHeight*(tan(omega1+theta1)+tan(omega1-theta1))/2;
328                     minRadius = kHeight*tan(omega1);
329                     eccentricity = sqrt(1-(minRadius*minRadius)/(maxRadius*maxRadius));
330                 
331
332
333                     offset = kHeight*(tan(omega1+theta1)-tan(omega1-theta1))/2;
334                     
335                     //printf("phi:%f theta:%f omega:%f \n", phi,theta,omega);
336
337                     //printf("offset:%f cx:%f cy:%f \n", offset,cx,cy);
338                      
339                     Float_t cxn = cx + offset * sin(phi);
340                     Float_t cyn = cy + offset * cos(phi);
341
342                     //printf("cxn:%f cyn:%f\n", cxn, cyn);
343
344                     for (Int_t cycle=0;cycle<maxcycle;cycle++)
345                       { 
346                         //printf("Entering point cycle");
347                         if(type==0)
348                           {
349                             points=(AliRICHDigit*) pDigits->UncheckedAt(cycle);
350                             segmentation->GetPadC(points->PadX(), points->PadY(),x, y, z);
351                           }
352                         
353                         if(type==1)
354                           {
355                             cluster=(AliRICHRawCluster*) pClusters->UncheckedAt(cycle);
356                             x=cluster->fX;
357                             y=cluster->fY;
358                             q=cluster->fQ;
359                           }
360                         
361                         if(type ==0 || q > 100)
362                           
363                           {
364
365                             x=x-cxn;
366                             y=y-cyn;
367                             radius=TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
368
369                             phi_relative = asin(y/radius);
370                             phi_relative = TMath::Abs(phi_relative - phi);
371                             
372                             angularadius = maxRadius*sqrt((1-eccentricity*eccentricity)/(1-eccentricity*eccentricity*cos(phi_relative)*cos(phi_relative)));
373                         
374                             //printf("omega:%f min:%f rad:%f max:%f\n",omega, angularadius-sigma,radius,angularadius+sigma);
375
376
377                             if((angularadius-sigma)<radius && (angularadius+sigma)>radius)
378                               {
379                                 printf("omega:%f min:%f rad:%f max:%f\n",omega, angularadius-sigma,radius,angularadius+sigma);
380
381                                 bintheta=theta*kDimensionTheta/kMaxTheta;
382                                 binphi=phi*kDimensionPhi/kMaxPhi;
383                                 binomega=omega*kDimensionOmega/(kMaxOmega-kMinOmega);
384                                 
385                                 if(Int_t(bintheta+0.5)==Int_t(bintheta))
386                                   inttheta=Int_t(bintheta);
387                                 else
388                                   inttheta=Int_t(bintheta+0.5);
389                                 
390                                 if(Int_t(binomega+0.5)==Int_t(binomega))
391                                   intomega=Int_t(binomega);
392                                 else
393                               intomega=Int_t(binomega+0.5);
394                                 
395                                 if(Int_t(binphi+0.5)==Int_t(binphi))
396                                   intphi=Int_t(binphi);
397                                 else
398                                   intphi=Int_t(binphi+0.5);
399                                 
400                                 //printf("Point added at %d %d %d\n",inttheta,intphi,intomega);
401                                 
402                                 //if(type==0)
403                                 point[inttheta][intphi][intomega]+=1;
404                                 //if(type==1)
405                                 //point[inttheta][intphi][intomega]+=(Int_t)(q);
406                                 
407                                 //printf("Omega stored:%d\n",intomega);
408                                 Points->Fill(inttheta,intphi,intomega,(float) 1);
409                                 ThetaPhi->Fill(inttheta,intphi,(float) 1);
410                                 OmegaTheta->Fill(inttheta,intomega,(float) 1);
411                                 OmegaPhi->Fill(intphi,intomega,(float) 1);
412                                 //printf("Filling at %d %d %d\n",Int_t(theta*kDimensionTheta/kMaxTheta),Int_t(phi*kDimensionPhi/kMaxPhi),Int_t(omega*kDimensionOmega/kMaxOmega));
413                               }
414                                 //if(omega<kMaxOmega)point[Int_t(theta)][Int_t(phi)][Int_t(omega)]+=1;
415                           }
416                       }
417                   }
418               }
419             
420           }
421                 //printf("Used %d digits for theta %3.1f\n",counter1, theta*180/kPi);
422         
423         printf("\n");
424
425         meanradius=meanradius/counter;
426         //printf("Mean radius:%f, counter:%d\n",meanradius,counter);
427         rechit[5]=meanradius;
428         //printf("Used %d digits\n",counter1);
429         //printf("\n");
430
431         if(nev<2)
432           {
433             if(nev==0)
434               {
435                 fc1->cd(1);
436                 Points->Draw("colz");
437                 fc1->cd(2);
438                 ThetaPhi->Draw("colz");
439                 fc1->cd(3);
440                 OmegaTheta->Draw("colz");
441                 fc1->cd(4);
442                 OmegaPhi->Draw("colz");
443                 fc3->cd();
444                 DigitsXY->Draw("colz");
445               }
446             else
447               {
448                 //fc1->cd(1);
449                 //Points->Draw("same");
450                 //fc1->cd(2);
451                 //ThetaPhi->Draw("same");
452                 //fc1->cd(3);
453                 //OmegaTheta->Draw("same");
454                 //fc1->cd(4);
455                 //OmegaPhi->Draw("same");       
456               }
457           }
458         
459     
460         //SPOT execute twice
461         for(Int_t s=0;s<kSpot;s++)
462           {
463             printf("     Applying Spot algorithm, pass %d\n", s);
464         
465         //buffer copy
466         for(i=0;i<=kDimensionTheta;i++)
467           {
468             for(j=0;j<=kDimensionPhi;j++)
469               {
470                 for(k=0;k<=kDimensionOmega;k++)
471                   {
472                     point1[i][j][k]=point[i][j][k];     
473                   }
474               }
475           }
476
477         //SPOT algorithm                        
478         for(i=1;i<kDimensionTheta-1;i++)
479           {
480             for(j=1;j<kDimensionPhi-1;j++)
481               {
482                 for(k=1;k<kDimensionOmega-1;k++)
483                   {
484                     if((point[i][k][j]>point[i-1][k][j])&&(point[i][k][j]>point[i+1][k][j])&&
485                        (point[i][k][j]>point[i][k-1][j])&&(point[i][k][j]>point[i][k+1][j])&&
486                        (point[i][k][j]>point[i][k][j-1])&&(point[i][k][j]>point[i][k][j+1]))
487                       {
488                         //cout<<"SPOT"<<endl;
489                         //Execute SPOT on point                                                                                         
490                         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]));    
491                         point1[i-1][k][j]=Int_t(SPOTp*point[i-1][k][j]);
492                         point1[i+1][k][j]=Int_t(SPOTp*point[i+1][k][j]);
493                         point1[i][k-1][j]=Int_t(SPOTp*point[i][k-1][j]);
494                         point1[i][k+1][j]=Int_t(SPOTp*point[i][k+1][j]);
495                         point1[i][k][j-1]=Int_t(SPOTp*point[i][k][j-1]);
496                         point1[i][k][j+1]=Int_t(SPOTp*point[i][k][j+1]);
497                       }
498                   }
499               }
500           }
501         
502         //copy from buffer copy
503         counter1=0;
504         for(i=1;i<kDimensionTheta;i++)
505           {
506             for(j=1;j<kDimensionPhi;j++)
507               {
508                 for(k=1;k<kDimensionOmega;k++)
509                   {
510                     point[i][j][k]=point1[i][j][k];
511                     if(nev<20)
512                       {
513                         if(s==kSpot-1)
514                           {
515                             if(point1[i][j][k] != 0)
516                               {
517                                 SpotPoints->Fill(i,j,k,(float) point1[i][j][k]);
518                                 //printf("Random number %f\n",random->Rndm(2));
519                                 //if(random->Rndm() < .2)
520                                   //{
521                                 SpotThetaPhi->Fill(i,j,(float) point1[i][j][k]);
522                                 SpotOmegaTheta->Fill(i,k,(float) point1[i][j][k]);
523                                 SpotOmegaPhi->Fill(j,k,(float) point1[i][j][k]);
524                                     counter1++;
525                                   //}
526                                 //printf("Filling at %d %d %d value %f\n",i,j,k,(float) point1[i][j][k]);
527                               }
528                           }
529                       }
530                     //if(point1[i][j][k] != 0)
531                       //printf("Last transfer point: %d, point1, %d\n",point[i][j][k],point1[i][j][k]);
532                   }
533               }
534           }
535       }
536     
537     //printf("Filled %d cells\n",counter1);
538
539         if(nev<2)
540           {
541             if(nev==0)
542               {
543                 fc2->cd(1);
544                 SpotPoints->Draw("colz");
545                 fc2->cd(2);
546                 SpotThetaPhi->Draw("colz");
547                 fc2->cd(3);
548                 SpotOmegaTheta->Draw("colz");
549                 fc2->cd(4);
550                 SpotOmegaPhi->Draw("colz");
551               }
552             else
553               {
554                 //fc2->cd(1);
555                 //SpotPoints->Draw("same");
556                 //fc2->cd(2);
557                 //SpotThetaPhi->Draw("same");
558                 //fc2->cd(3);
559                 //SpotOmegaTheta->Draw("same");
560                 //fc2->cd(4);
561                 //SpotOmegaPhi->Draw("same");   
562               }
563           }
564     
565     
566         //Identification is equivalent to maximum determination
567         max=0;maxi=0;maxj=0;maxk=0;
568         
569         printf("     Proceeding to identification");
570         
571         for(i=0;i<kDimensionTheta;i++)
572           for(j=0;j<kDimensionPhi;j++)
573             for(k=0;k<kDimensionOmega;k++)
574               if(point[i][j][k]>max)
575                 {
576                   //cout<<"maxi="<<i*90/dimension<<" maxj="<<j*90/dimension<<" maxk="<<k*kMaxOmega/dimension*180/kPi<<" max="<<max<<endl;
577                 maxi=i;maxj=j;maxk=k;
578                 max=point[i][j][k];
579                 printf(".");
580                 //printf("Max Omega %d, Max Theta %d, Max Phi %d (%d counts)\n",maxk,maxi,maxj,max);
581                 }
582         printf("\n");
583  
584         Float_t FinalOmega; 
585         Float_t FinalTheta;
586         Float_t FinalPhi;
587
588    
589         if(max>activ_tresh)
590           {
591             FinalOmega = maxk*(kMaxOmega-kMinOmega)/kDimensionOmega; 
592             FinalTheta = maxi*kMaxTheta/kDimensionTheta;
593             FinalPhi = maxj*kMaxPhi/kDimensionPhi;
594             
595             FinalOmega += kMinOmega;
596           }
597         else
598           {
599             FinalOmega = 0; 
600             FinalTheta = 0;
601             FinalPhi = 0;
602             
603             printf("     Ambiguous data!\n");
604           }
605
606             
607
608         //printf("Detected angle for height %3.1f and for center %3.1f %3.1f:%f\n",h,cx,cy,maxk*kPi/(kDimensionTheta*4));
609         printf("     Indentified angles: cerenkov - %f, theta - %3.1f, phi - %3.1f (%f activation)\n", FinalOmega, FinalTheta*180/kPi, FinalPhi*180/kPi, max);
610         //printf("Detected angle for height %3.1f and for center %3.1f %3.1f:%f\n",kHeight,cx,cy,maxk);
611
612         AngleAct->Fill(FinalOmega, (float) max);
613         Activation->Fill(max, (float) 1);
614
615         //fscanf(omegas,"%f",&realomega);
616         //fscanf(thetas,"%f",&realtheta);
617         //printf("Real Omega: %f",realomega);                   
618         //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;               
619     
620         //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));
621     
622     /*for(j=0;j<np;j++)
623       pointpp(maxj*90/kDimensionTheta,maxi*90/kDimensionPhi,maxk*kMaxOmega/kDimensionOmega*180/kPi,cx,cy);//Generates a point on the elipse*/               
624
625
626     //Start filling rec. hits
627     
628     rechit[0] = FinalTheta;
629     rechit[1] = 90*kPi/180 + FinalPhi;
630     rechit[2] = FinalOmega;
631     rechit[3] = cx;
632     rechit[4] = cy;
633
634     //CreatePoints(FinalTheta, 270*kPi/180 + FinalPhi, FinalOmega, kHeight);
635        
636     //printf ("track %d, theta %f, phi %f, omega %f\n\n\n",track,rechit[0],rechit[1],rechit[2]);
637     //printf("rechit %f %f %f %f %f\n",rechit[0],rechit[1],rechit[2],rechit[3],rechit[4]);
638     //printf("Chamber:%d",nch);
639
640
641       }
642
643     else   //if no cerenkovs
644       
645       {
646         
647         rechit[0] = 0;
648         rechit[1] = 0;
649         rechit[2] = 0;
650         rechit[3] = 0;
651         rechit[4] = 0;
652         originalOmega = 0; 
653         originalTheta = 0;
654         originalPhi =0;
655       }
656
657     
658     // fill rechits
659     pRICH->AddRecHit3D(nch-1,rechit,originalOmega, originalTheta, originalPhi);
660     printf("track %d, theta r:%f o:%f, phi r:%f o:%f, omega r:%f o:%f cx:%f cy%f\n\n\n", track, rechit[0], originalTheta, rechit[1], originalPhi, rechit[2], originalOmega, cx, cy);
661
662   }
663
664   if(type==1)  //reco from clusters
665     {
666       pRICH->ResetRawClusters();
667       //Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
668       //gAlice->TreeR()->GetEvent(track);
669       //printf("Going to branch %d\n",track);
670       //gAlice->GetEvent(nev);
671     }
672
673         
674   //printf("\n\n\n\n");
675   gAlice->TreeR()->Fill();
676   TClonesArray *fRec;
677   for (i=0;i<kNCH;i++) {
678     fRec=pRICH->RecHitsAddress3D(i);
679     int ndig=fRec->GetEntriesFast();
680     printf ("Chamber %d, rings %d\n",i+1,ndig);
681   }
682
683   fc4->cd(1);
684   AngleAct->Draw();
685   fc4->cd(2);
686   Activation->Draw();
687
688   pRICH->ResetRecHits3D();
689
690   free_i3tensor(point,0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
691   free_i3tensor(point1,0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
692   
693
694 }
695
696 Int_t ***AliRICHDetectV1::i3tensor(long nrl, long nrh, long ncl, long nch, long ndl, long ndh)
697 // allocate a Int_t 3tensor with range t[nrl..nrh][ncl..nch][ndl..ndh] 
698 {
699   long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
700   Int_t ***t;
701   
702   int NR_END=1; 
703
704   // allocate pointers to pointers to rows 
705   t=(Int_t ***) malloc((size_t)((nrow+NR_END)*sizeof(Int_t**)));
706   if (!t) printf("allocation failure 1 in f3tensor()");
707   t += NR_END;
708   t -= nrl;
709   
710   // allocate pointers to rows and set pointers to them 
711   t[nrl]=(Int_t **) malloc((size_t)((nrow*ncol+NR_END)*sizeof(Int_t*)));
712   if (!t[nrl]) printf("allocation failure 2 in f3tensor()");
713   t[nrl] += NR_END;
714   t[nrl] -= ncl;
715   
716   // allocate rows and set pointers to them 
717   t[nrl][ncl]=(Int_t *) malloc((size_t)((nrow*ncol*ndep+NR_END)*sizeof(Int_t)));
718   if (!t[nrl][ncl]) printf("allocation failure 3 in f3tensor()");
719   t[nrl][ncl] += NR_END;
720   t[nrl][ncl] -= ndl;
721   
722   for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
723   for(i=nrl+1;i<=nrh;i++) {
724     t[i]=t[i-1]+ncol;
725     t[i][ncl]=t[i-1][ncl]+ncol*ndep;
726     for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
727   }
728   
729   // return pointer to array of pointers to rows 
730   return t;
731 }
732
733
734 void AliRICHDetectV1::free_i3tensor(int ***t, long nrl, long nrh, long ncl, long nch,long ndl, long ndh)
735 // free a Int_t f3tensor allocated by i3tensor()
736 {
737   int NR_END=1; 
738
739   free((char*) (t[nrl][ncl]+ndl-NR_END));
740   free((char*) (t[nrl]+ncl-NR_END));
741   free((char*) (t+nrl-NR_END));
742 }