]>
Commit | Line | Data |
---|---|---|
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 | |
41 | ClassImp(AliRICHDetect) | |
42 | //___________________________________________ | |
43 | AliRICHDetect::AliRICHDetect() : TObject() | |
44 | { | |
237c933d | 45 | |
46 | // Default constructor | |
47 | ||
2685bf00 | 48 | fc1 = 0; |
49 | fc2 = 0; | |
50 | fc3 = 0; | |
51 | ||
c1076715 | 52 | } |
53 | ||
54 | //___________________________________________ | |
55 | AliRICHDetect::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 | //___________________________________________ | |
76 | AliRICHDetect::~AliRICHDetect() | |
77 | { | |
c1076715 | 78 | |
ac6e04fc | 79 | // Destructor |
80 | ||
c1076715 | 81 | } |
82 | ||
83 | ||
9c77ea58 | 84 | void 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 | 695 | Float_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 | 711 | Int_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 | ||
748 | void 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 | ||
759 | Float_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 | ||
786 | Float_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 | //________________________________________________________________________________ | |
817 | void 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 | |
880 | 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) | |
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 | } |