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