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