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