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