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