1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
20 #include "AliRICHDetectV1.h"
22 #include "AliRICHPoints.h"
23 #include "AliRICHDetect.h"
24 #include "AliRICHDigit.h"
25 #include "AliRICHRawCluster.h"
26 #include "AliRICHSegmentationV0.h"
28 #include "TParticle.h"
39 ClassImp(AliRICHDetectV1)
42 //___________________________________________
43 AliRICHDetectV1::AliRICHDetectV1() : AliRICHDetect()
46 // Default constructor
54 //___________________________________________
55 AliRICHDetectV1::AliRICHDetectV1(const char *name, const char *title)
56 :AliRICHDetect(name,title)
59 TStyle *mystyle=new TStyle("Plain","mystyle");
60 mystyle->SetPalette(1,0);
64 fc1= new TCanvas("c1","Reconstructed points",50,50,300,350);
66 fc2= new TCanvas("c2","Reconstructed points after SPOT",370,50,300,350);
68 fc3= new TCanvas("c3","Used Digits",690,50,300,350);
69 fc4= new TCanvas("c4","Mesh activation data",50,430,600,350);
75 //___________________________________________
76 AliRICHDetectV1::~AliRICHDetectV1()
84 void AliRICHDetectV1::Detect(Int_t nev, Int_t type)
88 // Detection algorithm
91 //printf("Detection started!\n");
92 Float_t omega,omega1,theta1,x,y,q=0,z,cx,cy,max,radius=0,meanradius=0;
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;
101 AliRICH *pRICH = (AliRICH*)gAlice->GetDetector("RICH");
102 AliRICHSegmentationV0* segmentation;
103 AliRICHChamber* iChamber;
104 AliRICHGeometry* geometry;
106 iChamber = &(pRICH->Chamber(0));
107 segmentation=(AliRICHSegmentationV0*) iChamber->GetSegmentationModel();
108 geometry=iChamber->GetGeometryModel();
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)
114 const Float_t kPi=TMath::Pi();
116 const Float_t kHeight=geometry->GetRadiatorToPads(); //Distance from Radiator to Pads in centimeters
117 //printf("Distance to Pads:%f\n",kHeight);
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
122 const Int_t kDimensionTheta=2; //Matrix dimension for angle Detection
123 const Int_t kDimensionPhi=2;
124 const Int_t kDimensionOmega=50;
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;
134 const Float_t sigma=0.5; //half thickness of fiducial band in cm
136 Float_t rechit[6]; //Reconstructed point data
138 Int_t ***point = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
139 Int_t ***point1 = i3tensor(0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
141 steptheta=(kMaxTheta-kMinTheta)/kDimensionTheta;
142 stepphi=(kMaxPhi-kMinPhi)/kDimensionPhi;
143 stepomega=(kMaxOmega-kMinOmega)/kDimensionOmega;
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");
180 Int_t ntracks = (Int_t)pRICH->TreeH()->GetEntries();
182 Float_t trackglob[3];
185 //printf("Area de uma elipse com teta 0 e Omega 45:%f",Area(0,45));
189 for (track=0; track<ntracks;track++) {
191 pRICH->TreeH()->GetEvent(track);
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();
197 AliRICHhit *mHit = 0;
200 Int_t counter=0, counter1=0;
202 for(i=0;i<kDimensionTheta;i++)
204 for(j=0;j<kDimensionPhi;j++)
206 for(k=0;k<kDimensionOmega;k++)
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)));
218 Int_t ncerenkovs = pRICH->Cerenkovs()->GetEntriesFast();
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)
232 originalOmega += cherenkov;
233 //printf("%f\n",cherenkov);
237 originalOmega = originalOmega/counter;
239 //printf("Cerenkovs : %d\n",counter);
241 mHit = (AliRICHhit*) pHits->UncheckedAt(0);
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();
250 printf("\n--------------------------------------\n");
251 printf("Chamber %d, track %d\n", nch, track);
252 printf("Original omega: %f\n",originalOmega);
254 iChamber = &(pRICH->Chamber(nch-1));
256 printf("Nch:%d x:%f y:%f\n",nch,trackglob[0],trackglob[2]);
258 iChamber->GlobaltoLocal(trackglob,trackloc);
260 iChamber->LocaltoGlobal(trackloc,trackglob);
266 printf("cy:%f ", cy);
268 if(counter != 0) //if there are cerenkovs
271 AliRICHDigit *points = 0;
272 TClonesArray *pDigits = pRICH->DigitsAddress(nch-1);
274 AliRICHRawCluster *cluster =0;
275 TClonesArray *pClusters = pRICH->RawClustAddress(nch-1);
279 //digitize from digits
283 gAlice->TreeD()->GetEvent(0);
284 Int_t ndigits = pDigits->GetEntriesFast();
286 //printf("Got %d digits\n",ndigits);
289 //digitize from clusters
293 Int_t nent=(Int_t)gAlice->TreeR()->GetEntries();
294 gAlice->TreeR()->GetEvent(nent-1);
295 Int_t nclusters = pClusters->GetEntriesFast();
297 //printf("Got %d clusters\n",nclusters);
304 printf("Starting calculations\n");
305 printf(" Start Finish\n");
306 printf("Progress: ");
309 for(Float_t theta=0;theta<kMaxTheta;theta+=steptheta)
312 for(Float_t phi=0;phi<=kMaxPhi;phi+=stepphi)
314 for(omega=kMinOmega;omega<=kMaxOmega;omega+=stepomega)
316 //printf("Entering angle cycle\n");
317 omega1=SnellAngle(omega);
318 theta1=SnellAngle(theta);
320 maxRadius = kHeight*(tan(omega1+theta1)+tan(omega1-theta1))/2;
321 minRadius = kHeight*tan(omega1);
322 eccentricity = sqrt(1-(minRadius*minRadius)/(maxRadius*maxRadius));
326 offset = kHeight*(tan(omega1+theta1)-tan(omega1-theta1))/2;
328 //printf("phi:%f theta:%f omega:%f \n", phi,theta,omega);
330 //printf("offset:%f cx:%f cy:%f \n", offset,cx,cy);
332 Float_t cxn = cx + offset * sin(phi);
333 Float_t cyn = cy + offset * cos(phi);
335 //printf("cxn:%f cyn:%f\n", cxn, cyn);
337 for (Int_t cycle=0;cycle<maxcycle;cycle++)
339 //printf("Entering point cycle");
342 points=(AliRICHDigit*) pDigits->UncheckedAt(cycle);
343 segmentation->GetPadC(points->PadX(), points->PadY(),x, y, z);
348 cluster=(AliRICHRawCluster*) pClusters->UncheckedAt(cycle);
354 if(type ==0 || q > 100)
360 radius=TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
362 phi_relative = asin(y/radius);
363 phi_relative = TMath::Abs(phi_relative - phi);
365 angularadius = maxRadius*sqrt((1-eccentricity*eccentricity)/(1-eccentricity*eccentricity*cos(phi_relative)*cos(phi_relative)));
367 //printf("omega:%f min:%f rad:%f max:%f\n",omega, angularadius-sigma,radius,angularadius+sigma);
370 if((angularadius-sigma)<radius && (angularadius+sigma)>radius)
372 printf("omega:%f min:%f rad:%f max:%f\n",omega, angularadius-sigma,radius,angularadius+sigma);
374 bintheta=theta*kDimensionTheta/kMaxTheta;
375 binphi=phi*kDimensionPhi/kMaxPhi;
376 binomega=omega*kDimensionOmega/(kMaxOmega-kMinOmega);
378 if(Int_t(bintheta+0.5)==Int_t(bintheta))
379 inttheta=Int_t(bintheta);
381 inttheta=Int_t(bintheta+0.5);
383 if(Int_t(binomega+0.5)==Int_t(binomega))
384 intomega=Int_t(binomega);
386 intomega=Int_t(binomega+0.5);
388 if(Int_t(binphi+0.5)==Int_t(binphi))
389 intphi=Int_t(binphi);
391 intphi=Int_t(binphi+0.5);
393 //printf("Point added at %d %d %d\n",inttheta,intphi,intomega);
396 point[inttheta][intphi][intomega]+=1;
398 //point[inttheta][intphi][intomega]+=(Int_t)(q);
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));
407 //if(omega<kMaxOmega)point[Int_t(theta)][Int_t(phi)][Int_t(omega)]+=1;
414 //printf("Used %d digits for theta %3.1f\n",counter1, theta*180/kPi);
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);
429 Points->Draw("colz");
431 ThetaPhi->Draw("colz");
433 OmegaTheta->Draw("colz");
435 OmegaPhi->Draw("colz");
437 DigitsXY->Draw("colz");
442 //Points->Draw("same");
444 //ThetaPhi->Draw("same");
446 //OmegaTheta->Draw("same");
448 //OmegaPhi->Draw("same");
454 for(Int_t s=0;s<kSpot;s++)
456 printf(" Applying Spot algorithm, pass %d\n", s);
459 for(i=0;i<=kDimensionTheta;i++)
461 for(j=0;j<=kDimensionPhi;j++)
463 for(k=0;k<=kDimensionOmega;k++)
465 point1[i][j][k]=point[i][j][k];
471 for(i=1;i<kDimensionTheta-1;i++)
473 for(j=1;j<kDimensionPhi-1;j++)
475 for(k=1;k<kDimensionOmega-1;k++)
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]))
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]);
495 //copy from buffer copy
497 for(i=1;i<kDimensionTheta;i++)
499 for(j=1;j<kDimensionPhi;j++)
501 for(k=1;k<kDimensionOmega;k++)
503 point[i][j][k]=point1[i][j][k];
508 if(point1[i][j][k] != 0)
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)
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]);
519 //printf("Filling at %d %d %d value %f\n",i,j,k,(float) point1[i][j][k]);
523 //if(point1[i][j][k] != 0)
524 //printf("Last transfer point: %d, point1, %d\n",point[i][j][k],point1[i][j][k]);
530 //printf("Filled %d cells\n",counter1);
537 SpotPoints->Draw("colz");
539 SpotThetaPhi->Draw("colz");
541 SpotOmegaTheta->Draw("colz");
543 SpotOmegaPhi->Draw("colz");
548 //SpotPoints->Draw("same");
550 //SpotThetaPhi->Draw("same");
552 //SpotOmegaTheta->Draw("same");
554 //SpotOmegaPhi->Draw("same");
559 //Identification is equivalent to maximum determination
560 max=0;maxi=0;maxj=0;maxk=0;
562 printf(" Proceeding to identification");
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)
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;
573 //printf("Max Omega %d, Max Theta %d, Max Phi %d (%d counts)\n",maxk,maxi,maxj,max);
584 FinalOmega = maxk*(kMaxOmega-kMinOmega)/kDimensionOmega;
585 FinalTheta = maxi*kMaxTheta/kDimensionTheta;
586 FinalPhi = maxj*kMaxPhi/kDimensionPhi;
588 FinalOmega += kMinOmega;
596 printf(" Ambiguous data!\n");
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);
605 AngleAct->Fill(FinalOmega, (float) max);
606 Activation->Fill(max, (float) 1);
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;
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));
616 pointpp(maxj*90/kDimensionTheta,maxi*90/kDimensionPhi,maxk*kMaxOmega/kDimensionOmega*180/kPi,cx,cy);//Generates a point on the elipse*/
619 //Start filling rec. hits
621 rechit[0] = FinalTheta;
622 rechit[1] = 90*kPi/180 + FinalPhi;
623 rechit[2] = FinalOmega;
627 //CreatePoints(FinalTheta, 270*kPi/180 + FinalPhi, FinalOmega, kHeight);
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);
636 else //if no cerenkovs
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);
657 if(type==1) //reco from clusters
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);
667 //printf("\n\n\n\n");
668 gAlice->TreeR()->Fill();
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);
681 pRICH->ResetRecHits3D();
683 free_i3tensor(point,0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
684 free_i3tensor(point1,0,kDimensionTheta,0,kDimensionPhi,0,kDimensionOmega);
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]
692 long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
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()");
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()");
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;
715 for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
716 for(i=nrl+1;i<=nrh;i++) {
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;
722 // return pointer to array of pointers to rows
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()
730 nrh++;ndh++;nch++;//to remove warning
733 free((char*) (t[nrl][ncl]+ndl-NR_END));
734 free((char*) (t[nrl]+ncl-NR_END));
735 free((char*) (t+nrl-NR_END));