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 **************************************************************************/
17 #include "AliRICHClusterFinder.h"
20 #include "AliRICHMap.h"
21 #include "AliRICHSDigit.h"
22 #include "AliRICHDigit.h"
23 #include "AliRICHRawCluster.h"
24 #include "AliRICHParam.h"
33 static AliSegmentation *gSegmentation;
34 static AliRICHResponse* gResponse;
35 static Int_t gix[500];
36 static Int_t giy[500];
37 static Float_t gCharge[500];
39 static Bool_t gFirst=kTRUE;
40 static TMinuit *gMyMinuit ;
41 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t);
42 static Int_t gChargeTot;
44 ClassImp(AliRICHClusterFinder)
45 //__________________________________________________________________________________________________
46 AliRICHClusterFinder::AliRICHClusterFinder(AliRICH *pRICH)
48 Info("main ctor","Start.");
52 fSegmentation=Rich()->C(1)->GetSegmentationModel();
53 fResponse =Rich()->C(1)->GetResponseModel();
55 fDigits=0; fNdigits=0;
64 //__________________________________________________________________________________________________
65 void AliRICHClusterFinder::Decluster(AliRICHRawCluster *cluster)
66 {// Decluster algorithm
67 Info("Decluster","Start.");
68 if(cluster->fMultiplicity==1||cluster->fMultiplicity==2){//Nothing special for 1- and 2-clusters
69 if(fNPeaks != 0) {cluster->fNcluster[0]=fNPeaks; cluster->fNcluster[1]=0;}
70 AddRawCluster(*cluster);
72 }else if(cluster->fMultiplicity==3){// 3-cluster, check topology
73 Centered(cluster);// ok, cluster is centered and added in Centered()
74 }else{//4-and more-pad clusters
75 if(cluster->fMultiplicity<= fMaxClusterSize){
76 SplitByLocalMaxima(cluster);
80 //__________________________________________________________________________________________________
81 Bool_t AliRICHClusterFinder::Centered(AliRICHRawCluster *cluster)
82 {//Is the cluster centered?
85 dig= (AliRICHDigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
86 Int_t x[kMaxNeighbours], y[kMaxNeighbours], xN[kMaxNeighbours], yN[kMaxNeighbours];
87 Int_t nn=Rich()->Param()->Neighbours(dig->PadX(),dig->PadY(),x,y);
91 for (Int_t i=0; i<nn; i++){//neighbours loop
92 if(fHitMap->TestHit(x[i],y[i]) == kUsed){
99 if(nd==2){// cluster is centered !
101 cluster->fNcluster[0]=fNPeaks;
102 cluster->fNcluster[1]=0;
105 AddRawCluster(*cluster);
109 // Highest signal on an edge, split cluster into 2+1
110 // who is the neighbour ?
111 Int_t nind=fHitMap->GetHitIndex(xN[0], yN[0]);
112 Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
113 Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
115 AliRICHRawCluster cnew;
117 cnew.fNcluster[0]=-1;
118 cnew.fNcluster[1]=fNRawClusters;
120 cnew.fNcluster[0]=fNPeaks;
123 cnew.fMultiplicity=2;
124 cnew.fIndexMap[0]=cluster->fIndexMap[0];
125 cnew.fIndexMap[1]=cluster->fIndexMap[i1];
127 cnew.fClusterType=cnew.PhysicsContribution();
131 cluster->fMultiplicity=1;
132 cluster->fIndexMap[0]=cluster->fIndexMap[i2];
133 cluster->fIndexMap[1]=0;
134 cluster->fIndexMap[2]=0;
135 FillCluster(cluster);
137 cluster->fNcluster[0]=fNPeaks;
138 cluster->fNcluster[1]=0;
140 cluster->fClusterType=cluster->PhysicsContribution();
141 AddRawCluster(*cluster);
145 Warning("Centered","\n Completely screwed up %d !! \n",nd);
150 //__________________________________________________________________________________________________
151 void AliRICHClusterFinder::SplitByLocalMaxima(AliRICHRawCluster *c)
152 {// Split the cluster according to the number of maxima inside
153 Info("SplitbyLocalMaxima","Start.");
155 AliRICHDigit* dig[100], *digt;
156 Int_t ix[100], iy[100], q[100];
157 Float_t x[100], y[100];
158 Int_t i; // loops over digits
159 Int_t j; // loops over local maxima
160 Int_t mul=c->fMultiplicity;
161 // dump digit information into arrays
162 for (i=0; i<mul; i++){
163 dig[i]= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
164 ix[i]= dig[i]->PadX();
165 iy[i]= dig[i]->PadY();
166 q[i] = dig[i]->Signal();
167 AliRICHParam::Pad2Local(ix[i], iy[i], x[i], y[i]);
172 Int_t associatePeak[100];
175 Int_t xNei[kMaxNeighbours], yNei[kMaxNeighbours];
176 for (i=0; i<mul; i++) {
177 fSegmentation->Neighbours(ix[i], iy[i], &nn, xNei, yNei);
179 for (j=0; j<nn; j++) {
180 if (fHitMap->TestHit(xNei[j], yNei[j])==kEmpty) continue;
181 digt=(AliRICHDigit*) fHitMap->GetHit(xNei[j], yNei[j]);
182 if (digt->Signal() > q[i]) {
185 // handle special case of neighbouring pads with equal signal
186 } else if (digt->Signal() == q[i]) {
188 for (Int_t k=0; k<nLocal; k++) {
189 if (xNei[j]==ix[indLocal[k]] && yNei[j]==iy[indLocal[k]]){
195 } // loop over next neighbours
196 // Maxima should not be on the edge
201 } // loop over all digits
202 // If only one local maximum found but multiplicity is high take global maximum from the list of digits.
203 if (nLocal==1 && mul>5) {
205 for (i=0; i<mul; i++) {
215 if(nLocal==2) {// If number of local maxima is 2 try to fit a double gaussian
217 // Initialise global variables for fit
219 gSegmentation=fSegmentation;
220 gResponse =fResponse;
223 for (i=0; i<mul; i++) {
226 gCharge[i]=Float_t(q[i]);
228 if (gFirst) gMyMinuit = new TMinuit(5);
230 gMyMinuit->SetFCN(fcn);
231 gMyMinuit->mninit(5,10,7);
232 Double_t arglist[20];
234 // Set starting values
235 static Double_t vstart[5];
236 vstart[0]=x[indLocal[0]];
237 vstart[1]=y[indLocal[0]];
238 vstart[2]=x[indLocal[1]];
239 vstart[3]=y[indLocal[1]];
240 vstart[4]=Float_t(q[indLocal[0]])/Float_t(q[indLocal[0]]+q[indLocal[1]]);
241 // lower and upper limits
242 static Double_t lower[5], upper[5];
243 lower[0]=vstart[0]-AliRICHParam::PadSizeX()/2;
244 lower[1]=vstart[1]-AliRICHParam::PadSizeY()/2;
246 upper[0]=vstart[0]+AliRICHParam::PadSizeX()/2;
247 upper[1]=vstart[1]+AliRICHParam::PadSizeY()/2;
249 lower[2]=vstart[2]-AliRICHParam::PadSizeX()/2;
250 lower[3]=vstart[3]-AliRICHParam::PadSizeY()/2;
252 upper[2]=vstart[2]+AliRICHParam::PadSizeX()/2;
253 upper[3]=vstart[3]+AliRICHParam::PadSizeY()/2;
258 static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
261 gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],iErr);
262 gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],iErr);
263 gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],iErr);
264 gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],iErr);
265 gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],iErr);
266 // ready for minimisation
267 gMyMinuit->SetPrintLevel(-1);
268 gMyMinuit->mnexcm("SET OUT", arglist, 0, iErr);
272 gMyMinuit->mnexcm("SET NOGR", arglist, 0, iErr);
273 gMyMinuit->mnexcm("SIMPLEX", arglist, 0, iErr);
274 gMyMinuit->mnexcm("MIGRAD", arglist, 0, iErr);
275 gMyMinuit->mnexcm("EXIT" , arglist, 0, iErr);
277 Double_t xrec[2], yrec[2], qfrac;
279 Double_t epxz, b1, b2;
280 gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, iErr);
281 gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, iErr);
282 gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, iErr);
283 gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, iErr);
284 gMyMinuit->mnpout(4, chname, qfrac, epxz, b1, b2, iErr);
286 cout<<"xrex[0]="<<xrec[0]<<"yrec[0]="<<yrec[0]<<"xrec[1]="<<xrec[1]<<"yrec[1]="<<yrec[1]<<"qfrac="<<qfrac<<endl;
287 for (j=0; j<2; j++) { // One cluster for each maximum
288 AliRICHRawCluster cnew;
290 cnew.fNcluster[0]=-1;
291 cnew.fNcluster[1]=fNRawClusters;
293 cnew.fNcluster[0]=fNPeaks;
296 cnew.fMultiplicity=0;
297 cnew.fX=Float_t(xrec[j]);
298 cnew.fY=Float_t(yrec[j]);
300 cnew.fQ=Int_t(gChargeTot*qfrac);
302 cnew.fQ=Int_t(gChargeTot*(1-qfrac));
304 for (i=0; i<mul; i++) {
305 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
306 cnew.fContMap[cnew.fMultiplicity]=AliRICHParam::AssignChargeToPad(xrec[j],yrec[j],gix[i], giy[i]);
307 cnew.fMultiplicity++;
309 FillCluster(&cnew,0);
310 cnew.fClusterType=cnew.PhysicsContribution();
314 }//if 2 maximum in cluster
317 if (nLocal >2 || !fitted) {
318 // Check if enough local clusters have been found, if not add global maxima to the list
323 Warning("SplitByLocalMaxima","no local maximum found");
327 if (nPerMax > fNperMax) {
328 Int_t nGlob=mul/fNperMax-nLocal+1;
331 for (i=0; i<mul; i++) {
338 if (nnew==nGlob) break;
342 for (i=0; i<mul; i++) { // Associate hits to peaks
345 if (isLocal[i]) continue;
346 for (j=0; j<nLocal; j++) {
347 Int_t il=indLocal[j];
348 Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
349 +(y[i]-y[il])*(y[i]-y[il]));
351 if (d<dmin) { // Select nearest peak
355 } else if (d==dmin) { // If more than one take highest peak
364 // One cluster for each maximum
365 for (j=0; j<nLocal; j++) {
366 AliRICHRawCluster cnew;
368 cnew.fNcluster[0]=-1;
369 cnew.fNcluster[1]=fNRawClusters;
371 cnew.fNcluster[0]=fNPeaks;
374 cnew.fIndexMap[0]=c->fIndexMap[indLocal[j]];
375 cnew.fMultiplicity=1;
376 for (i=0; i<mul; i++) {
377 if (isLocal[i]) continue;
378 if (associatePeak[i]==j) {
379 cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
380 cnew.fMultiplicity++;
384 cnew.fClusterType=cnew.PhysicsContribution();
389 }//SplitByLocalMaxima(AliRICHRawCluster *c)
390 //__________________________________________________________________________________________________
391 void AliRICHClusterFinder::FillCluster(AliRICHRawCluster* c, Int_t flag)
392 {// Completes cluster information starting from list of digits
406 for (Int_t i=0; i<c->fMultiplicity; i++){
407 dig= (AliRICHDigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
410 Int_t q=dig->Signal();
411 if (dig->Physics() >= dig->Signal()) {
413 } else if (dig->Physics() == 0) {
415 } else c->fPhysicsMap[i]=1;
416 // peak signal and track list
418 if (q>c->fPeakSignal) {
420 c->fTracks[0]=dig->Hit();
421 c->fTracks[1]=dig->Track(0);
422 c->fTracks[2]=dig->Track(1);
425 if (c->fContMap[i] > fraction) {
426 fraction=c->fContMap[i];
428 c->fTracks[0]=dig->Hit();
429 c->fTracks[1]=dig->Track(0);
430 c->fTracks[2]=dig->Track(1);
434 AliRICHParam::Pad2Local(ix,iy,x,y);
440 } // loop over digits
445 c->fX=fSegmentation->GetAnod(c->fX);
447 // apply correction to the coordinate along the anode wire
450 AliRICHParam::Local2Pad(x,y,ix,iy);
451 AliRICHParam::Pad2Local(ix,iy,x,y);
452 Int_t isec=fSegmentation->Sector(ix,iy);
453 TF1* cogCorr = fSegmentation->CorrFunc(isec-1);
456 Float_t yOnPad=(c->fY-y)/fSegmentation->Dpy(isec);
457 c->fY=c->fY-cogCorr->Eval(yOnPad, 0, 0);
461 //__________________________________________________________________________________________________
462 void AliRICHClusterFinder::AddDigit2Cluster(Int_t i, Int_t j, AliRICHRawCluster &c)
463 {//Find clusters Add i,j as element of the cluster
464 Info("AddDigit2Cluster","Start with digit(%i,%i)",i,j);
466 Int_t idx = fHitMap->GetHitIndex(i,j);
467 AliRICHDigit* dig = (AliRICHDigit*) fHitMap->GetHit(i,j);
468 Int_t q=dig->Signal();
469 if(q>TMath::Abs(c.fPeakSignal)){
471 c.fTracks[0]=dig->Hit();
472 c.fTracks[1]=dig->Track(0);
473 c.fTracks[2]=dig->Track(1);
475 // Make sure that list of digits is ordered
476 Int_t mu=c.fMultiplicity;
479 if (dig->Physics() >= dig->Signal()) {
481 } else if (dig->Physics() == 0) {
483 } else c.fPhysicsMap[mu]=1;
486 for (Int_t ind=mu-1; ind>=0; ind--) {
487 Int_t ist=(c.fIndexMap)[ind];
488 Int_t ql=((AliRICHDigit*)fDigits->UncheckedAt(ist))->Signal();
490 c.fIndexMap[ind]=idx;
491 c.fIndexMap[ind+1]=ist;
499 if (c.fMultiplicity >= 50 ) {
500 Info("AddDigit2CLuster","multiplicity >50 %d \n",c.fMultiplicity);
503 Float_t x,y;// Prepare center of gravity calculation
504 AliRICHParam::Pad2Local(i,j,x,y);
505 c.fX+=q*x; c.fY+=q*y; c.fQ += q;
506 fHitMap->FlagHit(i,j);// Flag hit as taken
509 Int_t xList[4], yList[4]; // Now look recursively for all neighbours
510 for (Int_t iNei=0;iNei<Rich()->Param()->Neighbours(i,j,xList,yList);iNei++)
511 if(fHitMap->TestHit(xList[iNei],yList[iNei])==kUnused) AddDigit2Cluster(xList[iNei],yList[iNei],c);
512 }//AddDigit2Cluster()
513 //__________________________________________________________________________________________________
514 void AliRICHClusterFinder::FindRawClusters()
515 {//finds neighbours and fill the tree with raw clusters
516 Info("FindRawClusters","Start for Chamber %i.",fChamber);
520 fHitMap=new AliRICHMap(fDigits);
522 for(Int_t iDigN=0;iDigN<fNdigits;iDigN++){//digits loop
523 AliRICHDigit *dig=(AliRICHDigit*)fDigits->UncheckedAt(iDigN);
524 Int_t i=dig->PadX(); Int_t j=dig->PadY();
525 if(fHitMap->TestHit(i,j)==kUsed||fHitMap->TestHit(i,j)==kEmpty) continue;
528 c.fMultiplicity=0; c.fPeakSignal=dig->Signal();
529 c.fTracks[0]=dig->Hit();c.fTracks[1]=dig->Track(0);c.fTracks[2]=dig->Track(1);
530 c.fNcluster[0]=-1;// tag the beginning of cluster list in a raw cluster
532 AddDigit2Cluster(i,j,c);//form initial cluster
534 c.fX /= c.fQ; // center of gravity
535 //c.fX=fSegmentation->GetAnod(c.fX);
539 // Int_t ix,iy;// apply correction to the coordinate along the anode wire
540 // Float_t x=c.fX, y=c.fY;
541 // Rich()->Param()->Local2Pad(x,y,ix,iy);
542 // Rich()->Param()->Pad2Local(ix,iy,x,y);
543 // Int_t isec=fSegmentation->Sector(ix,iy);
544 // TF1* cogCorr=fSegmentation->CorrFunc(isec-1);
546 // Float_t yOnPad=(c.fY-y)/fSegmentation->Dpy(isec);
547 // c.fY=c.fY-cogCorr->Eval(yOnPad,0,0);
550 c.fNcluster[1]=fNRawClusters; c.fClusterType=c.PhysicsContribution();
556 c.fMultiplicity=0; for(int k=0;k<c.fMultiplicity;k++) c.fIndexMap[k]=0;//reset cluster object
559 Info("FindRawClusters","Stop.");
561 //__________________________________________________________________________________________________
562 void AliRICHClusterFinder::CalibrateCOG()
570 fSegmentation->GiveTestPoints(n, x, y);
571 for (i=0; i<n; i++) {
575 SinoidalFit(xtest, ytest, func);
576 if (func) fSegmentation->SetCorrFunc(i, new TF1(*func));
580 //__________________________________________________________________________________________________
581 void AliRICHClusterFinder::SinoidalFit(Float_t x, Float_t y, TF1 *func)
583 static Int_t count=0;
588 Float_t xg[kNs], yg[kNs], xrg[kNs], yrg[kNs];
589 Float_t xsig[kNs], ysig[kNs];
592 AliRICHParam::Local2Pad(x,y,ix,iy);
593 AliRICHParam::Pad2Local(ix,iy,x,y);
594 Int_t isec=fSegmentation->Sector(ix,iy);
596 Float_t xmin = x-Rich()->Param()->PadSizeX()/2;
597 Float_t ymin = y-Rich()->Param()->PadSizeY()/2;
599 // Integration Limits
600 Float_t dxI=Rich()->Param()->SigmaIntegration()*Rich()->Param()->ChargeSpreadX();
601 Float_t dyI=Rich()->Param()->SigmaIntegration()*Rich()->Param()->ChargeSpreadY();
611 Float_t dy=Rich()->Param()->PadSizeY()/(kNs-1);
613 for (i=0; i<kNs; i++) {// Pad Loop
616 fSegmentation->SigGenInit(x, yscan, 0);
618 for (fSegmentation->FirstPad(x, yscan,0, dxI, dyI);
619 fSegmentation->MorePads();
620 fSegmentation->NextPad())
622 qp=fResponse->IntXY(fSegmentation);
626 Int_t ixs=fSegmentation->Ix();
627 Int_t iys=fSegmentation->Iy();
629 AliRICHParam::Pad2Local(ixs,iys,xs,ys);
633 Float_t ycog=sum/qcheck;
634 yg[i]=(yscan-y)/fSegmentation->Dpy(isec);
635 yrg[i]=(ycog-y)/fSegmentation->Dpy(isec);
641 Float_t dx=fSegmentation->Dpx(isec)/(kNs-1);
643 for (i=0; i<kNs; i++) {// Pad Loop
646 fSegmentation->SigGenInit(xscan, y, 0);
648 for (fSegmentation->FirstPad(xscan, y, 0, dxI, dyI);
649 fSegmentation->MorePads();
650 fSegmentation->NextPad())
652 qp=fResponse->IntXY(fSegmentation);
656 Int_t ixs=fSegmentation->Ix();
657 Int_t iys=fSegmentation->Iy();
659 AliRICHParam::Pad2Local(ixs,iys,xs,ys);
663 Float_t xcog=sum/qcheck;
664 xcog=fSegmentation->GetAnod(xcog);
666 xg[i]=(xscan-x)/fSegmentation->Dpx(isec);
667 xrg[i]=(xcog-x)/fSegmentation->Dpx(isec);
671 // Creates a Root function based on function sinoid above and perform the fit
672 TGraph *graphyr= new TGraph(kNs,yrg,ysig);
673 Double_t sinoid(Double_t *x, Double_t *par);
674 new TF1("sinoidf",sinoid,0.5,0.5,5);
675 graphyr->Fit("sinoidf","Q");
676 func = (TF1*)graphyr->GetListOfFunctions()->At(0);
678 //__________________________________________________________________________________________________
679 Double_t sinoid(Double_t *x, Double_t *par)
682 Double_t arg = -2*TMath::Pi()*x[0];
683 Double_t fitval= par[0]*TMath::Sin(arg)+
684 par[1]*TMath::Sin(2*arg)+
685 par[2]*TMath::Sin(3*arg)+
686 par[3]*TMath::Sin(4*arg)+
687 par[4]*TMath::Sin(5*arg);
690 //__________________________________________________________________________________________________
691 Double_t DoubleGauss(Double_t *x, Double_t *par)
692 {//Double gaussian function
693 Double_t arg1 = (x[0]-par[1])/0.18;
694 Double_t arg2 = (x[0]-par[3])/0.18;
695 return par[0]*TMath::Exp(-arg1*arg1/2)+par[2]*TMath::Exp(-arg2*arg2/2);
697 //__________________________________________________________________________________________________
698 Float_t DiscrCharge(Int_t i,Double_t *par)
700 // par[0] x-position of first cluster
701 // par[1] y-position of first cluster
702 // par[2] x-position of second cluster
703 // par[3] y-position of second cluster
704 // par[4] charge fraction of first cluster
705 // 1-par[4] charge fraction of second cluster
710 for (Int_t jbin=0; jbin<gNbins; jbin++) {
714 gChargeTot=Int_t(qtot);
717 Float_t q1=AliRICHParam::AssignChargeToPad(par[0],par[1],gix[i],giy[i]);
719 Float_t q2=AliRICHParam::AssignChargeToPad(par[2],par[3],gix[i],giy[i]);
720 // cout<<"qtot="<<gChargeTot<<" q1="<<q1<<" q2="<<q2<<" px="<<gix[i]<<" py="<<giy[i]<<endl;
722 Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
724 }//DiscrCharge(Int_t i,Double_t *par)
725 //__________________________________________________________________________________________________
726 void fcn(Int_t &npar, Double_t */*gin*/, Double_t &f, Double_t *par, Int_t)
727 {// Minimisation function
735 for (i=0; i<gNbins; i++) {
736 Float_t q0=gCharge[i];
737 Float_t q1=DiscrCharge(i,par);
738 delta=(q0-q1)/TMath::Sqrt(q0);
743 // chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
746 //__________________________________________________________________________________________________
747 void AliRICHClusterFinder::Exec()
749 Info("Exec","Start.");
751 Rich()->GetLoader()->LoadDigits();
753 for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events loop
754 gAlice->GetRunLoader()->GetEvent(iEventN);
756 Rich()->GetLoader()->MakeTree("R"); Rich()->MakeBranch("R");
757 Rich()->ResetDigitsOld(); Rich()->ResetRawClusters();
759 Rich()->GetLoader()->TreeD()->GetEntry(0);
760 for(fChamber=1;fChamber<=kNCH;fChamber++){//chambers loop
761 fDigits=Rich()->DigitsOld(fChamber); fNdigits=fDigits->GetEntries();
767 Rich()->GetLoader()->TreeR()->Fill();
768 Rich()->GetLoader()->WriteRecPoints("OVERWRITE");
770 Rich()->GetLoader()->UnloadDigits(); Rich()->GetLoader()->UnloadRecPoints();
771 Rich()->ResetDigitsOld(); Rich()->ResetRawClusters();
772 Info("Exec","Stop.");
774 //__________________________________________________________________________________________________