cuts on Q out, side, long added
[u/mrichter/AliRoot.git] / PMD / PMD_dig_reco.C
1 #include "PMD_dig_reco.h"
2 #include <iostream.h> // define cout stream
3 #include <fstream.h>  // defines ofstream class
4 #include <stdlib.h>   // defines exit() functions
5 void PMD_dig_reco (Int_t evNumber=1) 
6 {
7 /////////////////////////////////////////////////////////////////////////
8 //   This macro is a small example of a ROOT macro
9 //   illustrating how to read the output of GALICE
10 //   and fill some histograms.
11 //   
12 //
13 //  Created  by : Baba Potukuchi.V.K.S. ( JU / CERN ) 17/10/2000.  
14 //  Modified by : Baba Potukuchi.V.K.S. ( JU / CERN )  6/02/2001.
15 //  e-mail : potukuchi.baba@cern.ch 
16 //
17 /////////////////////////////////////////////////////////////////////////
18 // =========== Variables for Reconstruction==============================
19   FILE *fp4;
20   PMD_dig_reco* PPP = new PMD_dig_reco;
21 // ======================================================================
22   Int_t padmult,cpvmult,padmult_gam,tot_trk,tot_trk_cpv;
23   Int_t padmult_cut, cpvmult_cut;
24   Int_t trk_Mult[100000];
25   Int_t trk_Mult_cpv[100000];
26   Float_t xPos,yPos,zPos,Edep,Vol1,Vol2,Vol2s,Vol3,Vol4,Vol5;
27   Float_t  xpad,ypad;
28   Int_t iVol4,iVol5 , ixpad,iypad;
29   Int_t nbytes = 0;
30   //  Float_t px,py,pz;
31   Float_t px1,py1,pz1;
32   Float_t mass;
33   Float_t edep_tot;
34   Int_t id_part;
35   Int_t gam_in;
36   Int_t jjj;
37   Int_t pip_in,pim_in;
38   TObjArray  *points;
39   long double r,eta,px,py,pz;
40   Float_t phi,Eparticle,theta1,theta;
41   FILE *fp1;
42   FILE *fp2;
43   
44   fp1=fopen("kine_vertex.out","w");
45   fp2=fopen("kine_pmdhit.out","w");
46   fp4=fopen("det1.out","w");
47
48 // Dynamical link some shared libs
49    if (gClassTable->GetID("AliRun") < 0) {
50       gROOT->LoadMacro("loadlibs.C");
51       loadlibs();
52    }
53       
54 // Connect the Root Galice file containing Geometry, Kine and Hits
55    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
56    if (!file) file = new TFile("galice.root");
57
58 // Get AliRun object from file or create it if not on file
59 //   if (!gAlice) {
60       gAlice = (AliRun*)file->Get("gAlice");
61       if (gAlice) printf("AliRun object found on file\n");
62       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
63 //   }
64    // Just for test
65 //Loop over the requested number of events
66   for (Int_t iEvt = 0; iEvt < evNumber; iEvt++) {
67
68   printf("evt no %d \n",iEvt); 
69    // End Just for test
70
71
72
73 // Import the Kine and Hits Trees for the event evNumber in the file
74    Int_t nparticles = gAlice->GetEvent(iEvt);
75
76     printf ("N of particles %d\n",nparticles);
77
78     if (nparticles <= 0) return;
79
80         Int_t nevt = iEvt + 1;
81
82 //         ----********---- vertex track search BEGIN
83
84 //Loop over the requested number of events
85 // Import the Kine and Hits Trees for the event evNumber in the file
86
87         TObjArray *Particles = gAlice->Particles();
88         Float_t phi;
89    Float_t pl,cutmax,etamin,etamax,pmom,smin,smax;
90         Float_t *pxyz;
91         TObjArray  *points;
92         TParticle *particle;
93         AliPoints  *pm;
94         //     for (Int_t it = 0; it < 100000; it++)trk_Mult[it]=0.; 
95
96      for (Int_t in_track = 0; in_track < nparticles; in_track++) {
97      //     for (Int_t in_track = 0; in_track < 10; in_track++) {
98          particle = (TParticle*)gAlice->Particle(in_track);
99          if (!particle) continue;
100          id_part  = particle->GetPdgCode();
101          //      mass = TDatabasePDG::Instance()->GetParticle(id_part)->Mass();
102          //          printf("_____________--------%d %f\n",id_part,mass);
103
104          px  = particle->Px();
105          py  = particle->Py();
106          pz  = particle->Pz();
107          if(px !=0 || py!=0 || pz!=0) {
108          if( (px!=0 || py!=0) && pz < 0 ){
109 //         if(id_part == 22)gam_in++;
110          //        printf("%d\n",gam_in);         
111           r     = TMath::Sqrt(px*px + py*py);
112          theta = TMath::ATan2(r,TMath::Abs(pz));
113          phi = TMath::ATan2(py,px);
114          phi = (phi*180.)/3.14;
115 //         if(theta < -1.5707961)printf( " THETA=%f\n",theta) ;
116          if(theta > -1.5707961)eta = -TMath::Log(TMath::Tan(0.5*theta));
117          theta1=theta*180./3.14159;    
118
119 //         if(TMath::Abs(eta) > 2.3 && TMath::Abs(eta) < 3.5){
120          if(id_part == 22)gam_in++;
121          if(id_part == 211)pip_in++;
122          if(id_part == -211)pim_in++;
123      fprintf(fp1," %d %d %f %f %f %f %f %f \n",nevt,id_part,px,py,pz,eta,theta1,phi);
124 //            }
125 //         printf("%d %d %d\n",gam_in,pip_in,pim_in);         
126 //          printf("%f %f %f %f %d\n",ptg,pxg,pyg,pzg,id_part);
127      }
128     }
129   }
130         
131 // 1999 file-- Get pointers to Alice detectors and Hits containers
132      //------********----- vertex track ENDs
133
134    Float_t x,y,z,mass,e;
135    Int_t nbytes = 0;
136    Int_t j,hit,ipart;
137    Int_t npmd;
138    Int_t mbytes = 0;
139    Int_t nhits,gam_in;
140    Int_t ipmd,id_part;
141    Int_t ypad1;
142    Int_t mpmd;
143    const Int_t Npad2 = 72;
144
145    Int_t dVol;
146    //   Int_t iEvt;
147
148    Int_t sector,plane;
149    TParticle *particle;
150    AliPMDhit  *pmdHit;
151
152 // Get pointers to Alice detectors and Hits containers
153    AliDetector *PMD  = gAlice->GetDetector("PMD");
154 //   TClonesArray *fParticles = gAlice->Particles();
155    if (PMD) TClonesArray *PMDhits   = PMD->Hits();
156    if (PMD) TObjArray *points    = PMD->Points();
157
158    TTree *TH = gAlice->TreeH();
159    //   TTree *TK = gAlice->TreeK();
160
161    Int_t ntracks    = TH->GetEntries();
162    printf("No.of Tracksin the TreeH %d \n", ntracks);
163
164 // ===========   Initializing the array ============================
165  for (Int_t iDet = 1; iDet <= 2 ; ++iDet) {
166       printf(" Initializing for Detector %d \n",iDet) ;
167    for (Int_t iBox = 1; iBox <= 30 ; ++iBox) {
168      for (Int_t iCol = 1; iCol <= 72 ; ++iCol) {
169         for (Int_t iRow = 1; iRow <= 72 ; ++iRow) {
170        if(iDet==1)(*PPP).module1 [iBox] [iCol] [iRow] = 0.0 ;
171        if(iDet==2)(*PPP).module2 [iBox] [iCol] [iRow] = 0.0 ;
172         }
173      }
174    }
175  }
176 //  Creating histograms .
177    TH1F *edp = new TH1F("edp","Energy Deposited in each cell" ,100,100.,100000.);
178    TH1F *Modnum = new TH1F("Modnum","Module number" ,100,1.,100.);
179    TH2F *xycells = new TH2F("xycells","XY-of Cells in an Unit Modlule " ,100,0.,80.,100,0.,80.);
180    TH2F *xyhits = new TH2F("xyhits","Hits of the cells " ,100,-150.,150.,100,-150.,150.);
181
182 //  ================================================================
183     
184 // =======================================================
185 // Start loop on tracks in the hits containers
186    for (Int_t track=0; track<ntracks;track++) {
187         gAlice->ResetHits();
188
189      //changes introduced ******
190
191     nbytes += TH->GetEvent(track);
192     //   mbytes += TK->GetEvent(track);
193     //    printf("%d ******** %d \n" ,track,nbytes);
194      // ======>Histogram PMD
195
196      if (PMD) {
197        npmd = PMDhits->GetEntriesFast();
198        //printf("npmd=%d,  for track=%d %d \n" ,npmd,track,nbytes);
199          for (ipmd = 0; ipmd < npmd; ipmd++) {
200          pmdHit = (AliPMDhit*) PMDhits->UncheckedAt(ipmd);
201          ipart = pmdHit->GetTrack();
202       
203       //      printf("ID part and hits on PMD --------%d %d \n" ,ipart,track);      
204 //  get kinematics of the particles
205 //       particle = (TParticle*)fParticles->UncheckedAt(ipart);
206           particle = gAlice->Particle(ipart);
207 //       if (!particle) continue; 
208          //
209          id_part  = particle->GetPdgCode();
210          mass = TDatabasePDG::Instance()->GetParticle(id_part)->Mass();
211          //      printf("ID part and hit on PMD %d %d \n" ,id_part,ipmd);      
212          px  = particle->Px();
213          py  = particle->Py();
214          pz  = particle->Pz();
215          Eparticle  = particle->Energy();
216          
217          r     = TMath::Sqrt(px*px + py*py);
218          theta = TMath::ATan2(r,TMath::Abs(pz));
219          //phi = TMath::ATan2(TMath::Abs(py),TMath::Abs(px));
220          phi = TMath::ATan2(py,px);
221          phi = (phi*180.)/3.14159;
222          eta = -TMath::Log(TMath::Tan(0.5*theta));
223          if (pz < 0) eta = -eta;
224          theta=theta*180./3.14159;
225          
226 //         if(TMath::Abs(eta) > 2.3 && TMath::Abs(eta) < 3.5){
227          //
228          xPos = pmdHit->X();
229          yPos = pmdHit->Y();
230          zPos = pmdHit->Z();
231          Vol1 = pmdHit->fVolume[0];
232          Vol2 = pmdHit->fVolume[1];
233          Vol3 = pmdHit->fVolume[2];
234          Vol4 = pmdHit->fVolume[3];
235          Vol5 = pmdHit->fVolume[4];
236                  
237          Vol2s = Vol2-1;
238          ypad = Vol2s/Npad2 + 1;
239          ypad1 = ypad-1;
240          xpad = Vol2-ypad1*Npad2;
241          Edep = pmdHit->fEnergy;
242          ixpad = xpad+0.2 ;
243          iypad = ypad+0.2 ;
244          iVol4 = Vol4+0.2 ;
245          iVol5 = Vol5+0.2 ;
246 //         module [Vol4] [Vol5] [xpad] [ypad] = Edep ;
247          if(Edep>0)printf(" VOL4 %d %d %d %d %f \n",iVol4,iVol5,ixpad,iypad,Edep) ; 
248          edp ->Fill(Edep) ;
249          Modnum ->Fill(iVol5) ;
250          xycells->Fill(ixpad,iypad);
251          xyhits->Fill(xPos,yPos);
252         
253        if(iVol4==1)(*PPP).module1 [iVol5] [ixpad] [iypad] = Edep  ;
254        if(iVol4==2)(*PPP).module2 [iVol5] [ixpad] [iypad] = Edep  ;
255           Float_t XXP = (*PPP).module2 [iVol5] [ixpad] [iypad] ;
256          if(Edep>0)printf(" VOL4 %d %d %d %d %f \n",Vol4,Vol5,xpad,ypad,XXP) ; 
257          // Write to the kine.out file
258          fprintf(fp2," %d %d %d %f %f %f %d %d %d %d %f %f %f %f %f \n",
259                 nevt,ipart,id_part,xPos,yPos,zPos,xpad,ypad,Vol5,Vol4,Edep,px,py,pz,Eparticle);
260
261 //       Closing Eta loop of the incident particles.  
262 //         }
263          
264        }
265      }
266    }
267   }   
268
269   pmdreco_(PPP);
270
271 // TCanvas *c1 = new TCanvas("c1","PMD Simulated digits(Edep) ",400,10,600,700);
272 // c1->Divide(1,2);
273 // c1->cd(1);
274 // gPad->SetFillColor(33);
275 // edp->SetFillColor(42);
276 // edp->Draw();
277 // c1->cd(2);
278 // gPad->SetFillColor(33);
279 // Modnum->SetFillColor(41);
280 // Modnum->Draw();
281 // c1->Print("baba.ps");
282 // TCanvas *c2 = new TCanvas("c2","PMD Simulated digits(Edep) ");
283 // gPad->SetFillColor(33);
284 // xycells->SetFillColor(42);
285 // xycells->Draw();
286 // c2->Print("baba1.ps");
287 // TCanvas *c3 = new TCanvas("c3","PMD Simulated digits(Edep) ");
288 // gPad->SetFillColor(33);
289 // xyhits->SetFillColor(42);
290 // xyhits->Draw();
291 // c3->Print("baba2.ps");
292 }
293
294 // ===========  RECONSTRUCTION ======================================
295 void pmdreco_(PMD_dig_reco* PPP) {
296     Float_t  sqrth;
297     sqrth = sqrt(3.) / 2.;
298     
299     printf(" I am in MAIN calling RanMAR & Intpts_ \n") ;
300     rmarin_(PPP);
301     intpts_(PPP);
302 //  Generate COORDINATES .
303         for (Int_t i__ = 1; i__ <= 72; ++i__) {
304             for (Int_t j = 1; j <= 72; ++j) {
305                 (*PPP).coord [1] [i__] [j] = i__ + j / 2.;
306                 (*PPP).coord [2] [i__] [j] = j * sqrth;
307             }
308         }
309 //-------------------------
310   for (Int_t iDet = 1; iDet <= 2 ; ++iDet) {
311        printf("idet= %d \n",iDet) ;
312     for (Int_t iBox = 1; iBox <= 30 ; ++iBox) {
313       for (Int_t iCol = 1; iCol <= 72 ; ++iCol) {
314          for (Int_t iRow = 1; iRow <= 72 ; ++iRow) {
315           Float_t PVAL = (*PPP).module1 [iBox] [iCol] [iRow];
316 //       printf("Det1 = %f %f %f %f %f \n",iDet, iBox, iCol, iRow,PVAL ) ;
317           if(iDet==1)(*PPP).det1 [iCol] [iRow] = (*PPP).module1 [iBox] [iCol] [iRow] ;
318           if(iDet==2)(*PPP).det2 [iCol] [iRow] = (*PPP).module2 [iBox] [iCol] [iRow] ;
319        if(PVAL>0.0)printf("Det1 = %d %d %d %d %f \n",iDet, iBox, iCol, iRow,PVAL ) ;
320 //       fprintf(fp4," Det1 = %d %d %d %d %f \n",iDet, iBox, iCol, iRow,VAL ) ;
321
322          }
323       }
324 // Ending Reading pads in the Box .
325     }
326 // Ending Reading pads for 1-CPV,2-PMD.
327  }
328 /*     analyse detector 1 */
329 /*     CPV clustering. Ordering data */
330        printf(" I am in MAIN Doing Clustering for CPV.  \n") ;
331        printf (" I am in MAIN calling ==order== for CPV(det=1) \n");
332        order_(1,PPP);
333     printf(" I am in MAIN  after order for CPV(det=1) \n") ;
334     printf(" I am in MAIN calling ==crclust== \n") ;
335         crclust_(1, PPP) ;
336     printf(" I am in MAIN after ==crclust==  & now call prints \n") ;
337         prints_(1,PPP) ;
338        printf(" I am in MAIN Doing Clustering for PMD.  \n") ;
339 printf (" I am in MAIN calling ==order== for PMD(det=2) \n");
340         order_(2,PPP);
341     printf(" I am in MAIN  after order for PMD(det=2) \n") ;
342     printf(" I am in MAIN calling ==crclust== \n") ;
343         crclust_(2, PPP) ;
344     printf(" I am in MAIN after ==crclust==  & now call prints \n") ;
345         prints_(2,PPP) ;
346         (*PPP).incr = 0;
347         for (Int_t i__ = 1; i__ <= 72; ++i__) {
348             for (j = 1; j <= 72; ++j) {
349                 if ((*PPP).infocl[1] [i__] [j]  == 1) {
350                     ++(*PPP).incr;
351                     (*PPP).inford[1] [(*PPP).incr] = (*PPP).infocl[2] [i__] [j] ;
352                     (*PPP).inford[2] [(*PPP).incr] = i__;
353                     (*PPP).inford[3] [(*PPP).incr] = j;
354                 }
355             }
356         }
357         for (i__ = 1; i__ <= 72; ++i__) {
358             for (j = 1; j <= 72; ++j) {
359                 if ((*PPP).infocl[1] [i__] [j]  == 2) {
360                     ++(*PPP).incr;
361                     (*PPP).inford[1] [(*PPP).incr] = (*PPP).infocl[2] [i__] [j];
362                     (*PPP).inford[2] [(*PPP).incr] = i__;
363                     (*PPP).inford[3] [(*PPP).incr] = j;
364                 }
365             }
366         }
367         arrange_(PPP);
368 /*     better clustering ( Gaussian fit ) */
369         refclust_(PPP);
370     return ;
371 } /* LOOP pmdreco */
372 // 
373  /* >>>>>>>>>>>>  */void rmarin_(PMD_dig_reco* PPP)
374 {
375     Int_t i__, j, k, l, m;
376     Float_t s, t;
377 /*         Initializing routine for RANMAR, must be called before */
378 /*         generating any pseudorandom numbers with RANMAR. The input */
379 /*         values should be in the ranges:  0<=IJ<=31328 */
380 /*                                          0<=KL<=30081 */
381 /* To get the standard values in Marsaglia's paper, IJ=1802, KL=9373 */
382     Int_t ij = 8775;
383     Int_t kl = 5926;
384     if (ij <= 0 || ij >= 31328) {
385         goto L900;
386     }
387     if (kl <= 0 || kl >= 30081) {
388         goto L900;
389     }
390  printf("I am in RanMAR \n");
391     i__ = ij / 177 % 177 + 2;
392     j = ij % 177 + 2;
393     k = kl / 169 % 178 + 1;
394     l = kl % 169;
395     printf(" ranmar initialized: %d %d %d %d %d %d \n",ij,kl,i__,j,k,l);
396     for (Int_t ii = 1; ii <= 97; ++ii) {
397         s = 0.;
398         t = .5;
399         for (Int_t jj = 1; jj <= 24; ++jj) {
400             m = (((i__ * j) % 179) * k) % 179;
401             i__ = j;
402             j = k;
403             k = m ;
404             l = l * 53 + 1 % 169;
405             if (l * m % 64 >= 32) {
406                 s += t;
407             }
408             t *= .5;
409         }
410         (*PPP).u[ii - 1] = s;
411     }
412     (*PPP).c__ = .021602869033813477;
413     (*PPP).cd  = .45623308420181274;
414     (*PPP).cm  = .99999982118606567;
415      printf(" Exiting from rmarin %f %f %f \n",(*PPP).c__,(*PPP).cd,(*PPP).cm) ;
416      return ;
417 L900:
418     printf( " +++++++ stop in rmarin: wrong params! ");
419 } /* rmarin_ */
420
421 /*>>>>>>>>>>>>>>>> */void order_(Int_t jdet,PMD_dig_reco* PPP)
422 {
423     Int_t nmx,i__1, i__2, i__3;
424
425     Float_t  adum;
426     Int_t idum, i__, j, k, l;
427     Float_t  d1[5185];
428     Float_t  iord1[5185];
429     Int_t itest, i1, i2, id ;
430
431     printf("I am in order for DET=%d \n",jdet);
432     nmx = 5184;
433     for (i1 = 1; i1 <= 72; ++i1) {
434         for (i2 = 1; i2 <= 72; ++i2) {
435             i__ = (i2 - 1) * 72 + i1;
436             d1[i__] = 0.0 ;
437            iord1 [i__] = 0;
438            if(jdet==1) d1[i__] = (*PPP).det1 [i1] [i2] ;
439            if(jdet==2) d1[i__] = (*PPP).det2 [i1] [i2] ;
440            iord1 [i__] = i__;
441         }
442     }
443     for (j = 2; j <= nmx; ++j) {
444         itest = 0;
445         adum = d1[j];
446         idum = iord1 [j];
447         i__2 = j - 1;
448         for (k = 1; k <= i__2; ++k) {
449             if (adum > d1[k ] && itest == 0) {
450                 itest = 1 ; 
451                 i__3 = k;
452                 for (l = j - 1; l >= i__3; --l) {
453                     d1[l+1] = d1[l];
454                     iord1 [l+1] = iord1 [l];
455                 }
456                 d1[k] = adum;
457                 iord1 [k] = idum;
458             }
459         }
460     }
461     i__1 = nmx;
462     for (i__ = 1; i__ <= i__1; ++i__) {
463         id = iord1 [i__];
464         i2 = id / 72 + 1;
465         i1 = id - (i2 - 1) * 72;
466         (*PPP).iorder[1] [i__] = i1;
467         (*PPP).iorder[2] [i__] = i2;
468      
469      printf("I am in order for DET=%d %d %d %d %d %d \n",jdet,i__,i1,i2,(*PPP).iorder[1] [i__],(*PPP).iorder[2] [i__]);
470     }
471 } /* order_ */
472
473 /* crclust_ */ void crclust_(Int_t ndet, PMD_dig_reco* PPP)
474  { 
475            Float_t d__ [73] [73] ;
476
477     Int_t neibx[7] =  { 0,1,0,-1,-1,0,1 };
478     Int_t neiby[7] =  { 0,0,1,1,0,-1,-1 };
479     
480     Int_t i__1;
481     Float_t sqrth ;
482
483     Int_t icld, idum, i__, j, k;
484     Int_t nb, id1, id2, jd1, jd2, icle, icn, itr, nmx ;
485
486    printf(" neibx %d %d %d %d %d %d \n",neibx[1],neibx[2],neibx[3],neibx[4],neibx[5],neibx[6]) ;
487    printf(" neiby %d %d %d %d %d %d \n",neiby[1],neiby[2],neiby[3],neiby[4],neiby[5],neiby[6]) ;
488   
489         for (j = 1; j <= 72; ++j) {
490             for (k = 1; k <= 72; ++k) {
491              d__ [j] [k] = 0.0 ;
492              if(ndet==1) d__ [j] [k] = (*PPP).det1 [j] [k] ;
493              if(ndet==2) d__ [j] [k] = (*PPP).det2 [j] [k] ;
494              
495             }
496         }
497     nmx = 5148;
498     sqrth=TMath::Sqrt(3.)/2. ;
499 /*     initialise infocl */
500     for (i__ = 1; i__ <= 2; ++i__) {
501         for (j = 1; j <= 72; ++j) {
502             for (k = 1; k <= 72; ++k) {
503                 (*PPP).infocl[i__] [j] [k] = 0;
504             }
505         }
506     }
507 /*     compute number of pads with nonzero count ( nmx1 ) */
508     i__1 = nmx;
509     for (i__ = 1; i__ <= i__1; ++i__) {
510         id1 = (*PPP).iorder [1] [i__];
511         id2 = (*PPP).iorder [2] [i__];
512         if (d__[id1] [id2] > 0.0000000000001) {
513             (*PPP).nmx1 = i__;
514           printf(" nmx1 -> %d %d %f %d \n",id1,id2, d__ [id1] [id2],(*PPP).nmx1 ) ;
515         }
516     }
517     printf(" nmx1= %d \n",(*PPP).nmx1) ;
518 /*     crude clustering. start with pad with largest count ( cluster 1 ) 
519 */
520     icle = 1;
521     idum = 0;
522     id1 = (*PPP).iorder [1] [1];
523     id2 = (*PPP).iorder [2] [1];
524                 (*PPP).infocl[1] [id1] [id2] = 1;
525                 (*PPP).infocl[2] [id1] [id2] = icle;
526           printf(" infocl  %d %d %d %d %d \n",id1,id2,icle,(*PPP).infocl[1] [id1] [id2],(*PPP).infocl[2] [id1] [id2] ) ;
527 /*     check if neighbouring count in nonzero. If yes, include in the */
528 /*     cluster */
529 /*     infocl(1,i1,i2) =1 for primary ( peak ) =2 for secondary ( */
530 /*     neighbour ) */
531 /*     infocl(2,i1,i2) stores cluster no. */
532     for (nb = 1; nb <= 6; ++nb) {
533         jd1 = id1 + neibx[nb];
534         jd2 = id2 + neiby[nb];
535           printf(" infocl  %d %d %d %d %d %f %d %d\n",nb,jd1,id1,jd2,id2,d__[jd1] [jd2],icle,(*PPP).infocl[1] [jd1] [jd2] ) ;
536         if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
537             if (d__[jd1] [jd2] > 0.0) {
538                 (*PPP).infocl[1] [jd1] [jd2] = 2;
539                 (*PPP).infocl[2] [jd1] [jd2] = icle;
540             }
541         }
542           printf(" infocl  %d %d %d %d %d %f %d %d\n",nb,jd1,id1,jd2,id2,d__[jd1] [jd2],icle,(*PPP).infocl[1] [jd1] [jd2] ) ;
543     }
544 /*     start checking with next largest count and so on */
545     i__1 = (*PPP).nmx1;
546     for (i__ = 2; i__ <= i__1; ++i__) {
547         itr = 0;
548         id1 = (*PPP).iorder[1] [i__];
549         id2 = (*PPP).iorder[2] [i__];
550        printf("iorder %d %d %d %d %d %d %d\n",i__,id1,id2,(*PPP).infocl [1] [id1] [id2],jd1,jd2,(*PPP).infocl [1] [jd1] [jd2] ) ; 
551 /*     nonzero infocl(1,id1,id2) means the pad is included in earlier 
552 */
553 /*     cluster ( itr = 1 ) Otherwise itr = 0. */
554         if ((*PPP).infocl [1] [id1] [id2] > 0) {
555             itr = 1;
556         }
557         if (itr ==  0) {
558             for (nb = 1; nb <= 6; ++nb) {
559                 jd1 = id1 + neibx[nb ];
560                 jd2 = id2 + neiby[nb ];
561 /*     check if a neighbouring pad belongs to previous cluster
562 . If yes, */
563 /*     icn = 1. Otherwise icn = 0. */
564                 icn = 0;
565                 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
566                     if ((*PPP).infocl [1] [jd1] [jd2] != 0) {
567                         icn = 1;
568                         icld = (*PPP).infocl [2] [jd1] [jd2] ;
569                     }
570                 }
571             }
572 /*     itr = 0 and icn = 0 means new cluster. Store the cluster in
573 fo. in */
574 /*     infocl. */
575             if (TMath::Abs(icn) <= 0.1) {
576                 ++icle;
577                 (*PPP).infocl[1] [id1] [id2] = 1;
578                 (*PPP).infocl[2] [id1] [id2] = icle;
579               printf("ID1 icn %d %d %d %d %d\n",icn,icle,id1,id2,(*PPP).infocl[2] [id1] [id2]);
580                 for (nb = 1; nb <= 6; ++nb) {
581 /*     neighbours of the peak put with the peak. */
582                     jd1 = id1 + neibx[nb ];
583                     jd2 = id2 + neiby[nb ];
584                     if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
585             if (d__[jd1] [jd2] != 0  && (*PPP).infocl [1] [jd1] [jd2] == 0) {
586                 (*PPP).infocl[1] [jd1] [jd2] = 2;
587                 (*PPP).infocl[2] [jd1] [jd2] = icle;
588               printf("JD1 icn %d %d %d %d %d\n",icn,icle,jd1,jd2,(*PPP).infocl[2] [jd1] [jd2]);
589                         }
590                     }
591                 }
592            } else {
593 /*     itr = 0 and icn = 1 means the pad belongs to cluster ic
594 ld */
595                 (*PPP).infocl[1] [id1] [id2] = 2;
596                 (*PPP).infocl[2] [id1] [id2] = icld;
597                 for (nb = 1; nb <= 6; ++nb) {
598 /*     and also the neighbours. */
599                     jd1 = id1 + neibx[nb ];
600                     jd2 = id2 + neiby[nb ];
601                     if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
602                         if (d__[jd1] [jd2]  != 0. && (*PPP).infocl[1] [jd1] [jd2] == 0) {
603                         (*PPP).infocl[1] [jd1] [jd2] = 2;
604                         (*PPP).infocl[2] [jd1] [jd2] = icld;
605                         }
606                     }
607                 }
608          } else {
609 /*     itr = 0, so not a new cluster */
610             icld = (*PPP).infocl[2] [id1] [id2] ;
611             for (nb = 1; nb <= 6; ++nb) {
612 /*     neighbours put with the previous cluster. */
613                 jd1 = id1 + neibx[nb ];
614                 jd2 = id2 + neiby[nb ];
615                 if (jd1 > 0 && jd1 <= 72 && jd2 > 0 && jd2 <= 72) {
616                 if (d__[jd1] [jd2]  != 0. && (*PPP).infocl[1] [jd1] [jd2] == 0) {
617                         (*PPP).infocl[1] [jd1] [jd2] = 2;
618                         (*PPP).infocl[2] [jd1] [jd2] = icld;
619                     }
620                 }
621             }
622         }
623     }
624 }
625     printf(" ndet=%d,nmx1= %d , icle= %d icld= %d \n",ndet,(*PPP).nmx1,icle,icld) ;
626
627 /* L2: */
628 /* L3: */
629     return ;
630 } /* crclust_ */
631
632 /* Subroutine */ void intpts_(PMD_dig_reco* PPP)
633 {
634     Int_t i__1, i__2, i__3;
635
636     Int_t ixpt;
637     Float_t x, y, dx, dy;
638     Int_t ix, iy;
639     Float_t yl;
640     Int_t idx, idy, ipt;
641     Float_t wtx, wty;
642
643     ixpt = 2;
644     idx = -1;
645     ipt = 0;
646     dx = .25;
647     i__1 = ixpt;
648     for (ix = -ixpt; ix <= i__1; ++ix) {
649         x = ix * dx;
650         idx = -idx;
651         i__2 = ix / ixpt ;
652         wtx = (3. - idx -  (TMath::Abs(i__2) ))/ 3. * dx;
653         yl = 1. / TMath::Sqrt(3.) * (1. - TMath::Abs(x));
654         dy = yl * 2. / 4.;
655         idy = -1;
656         i__2 = ixpt;
657         for (iy = -ixpt; iy <= i__2; ++iy) {
658             idy = -idy;
659             y = iy * dy;
660             i__3 = iy / ixpt ;
661             wty = (3. - idy - (TMath::Abs(i__3))) / 3. * dy;
662             ++ipt;
663             (*PPP).pts[1] [ipt] = x;
664             (*PPP).pts[2] [ipt] = y;
665             (*PPP).wts[ipt] = wtx * wty;
666         }
667     }
668      printf(" Exiting from intpts_ %f %f  \n", ipt,(*PPP).wts[ipt] );
669     return ;
670 } /* intpts_ */
671
672 /* prints_  */ void prints_(Int_t ndet, PMD_dig_reco* PPP)
673 {
674     Int_t i__, j;
675
676     if (ndet == 1) {
677        FILE *fp1;
678        FILE *fp2;
679        fp1=fopen("CPV_prime.out","w");
680        fp2=fopen("CPV_secondary.out","w");
681     }
682     if (ndet == 2) {
683        FILE *fp1;
684        FILE *fp2;
685        fp1=fopen("PMD_prime.out","w");
686        fp2=fopen("PMD_secondary.out","w");
687     }
688     for (i__ = 1; i__ <= 72; ++i__) {
689         for (j = 1; j <= 72; ++j) {
690          if( (*PPP).infocl [1] [i__] [j] == 1) {
691           fprintf(fp1," %f %f %d %d %d \n",(*PPP).coord [1] [i__] [j],(*PPP).coord [2] [i__] [j],i__,j,(*PPP).infocl [2] [i__] [j] ) ;
692          }
693          if( (*PPP).infocl [1] [i__] [j] == 2) {
694           fprintf(fp2," %f %f %d %d %d \n",(*PPP).coord [1] [i__] [j],(*PPP).coord [2] [i__] [j],i__,j,(*PPP).infocl [2] [i__] [j] ) ;
695          }
696         }
697     }
698     return ;
699 } /* prints_ */
700
701 /* Subroutine */void  arrange_(PMD_dig_reco* PPP)
702 {
703     Int_t i__1, i__2, i__3;
704
705     Int_t ihld1, ihld2, ihld3, i__, j, k, itest;
706
707     i__1 = (*PPP).incr;
708     for (i__ = 2; i__ <= i__1; ++i__) {
709         itest = 0;
710         ihld1 = (*PPP).inford [1] [i__];
711         ihld2 = (*PPP).inford [2] [i__];
712         ihld3 = (*PPP).inford [3] [i__];
713         i__2 = i__ - 1;
714         for (j = 1; j <= i__2; ++j) {
715             if (itest == 0 && (*PPP).inford [1] [j] > ihld1) {
716                 itest = 1;
717                 i__3 = j;
718                 for (k = i__ -1; k >= i__3; --k) {
719                  (*PPP).inford [1] [k+1] = (*PPP).inford [1] [k] ;
720                  (*PPP).inford [2] [k+1] = (*PPP).inford [2] [k] ;
721                  (*PPP).inford [3] [k+1] = (*PPP).inford [3] [k] ;
722                 }
723                  (*PPP).inford [1] [j] = ihld1 ;
724                  (*PPP).inford [2] [j] = ihld2 ;
725                  (*PPP).inford [3] [j] = ihld3 ;
726             }
727         }
728     }
729      printf(" Exiting from arrange \n") ;
730     return ;
731 } /* arrange_ */
732
733 /* refclust_ */void  refclust_(PMD_dig_reco* PPP)
734 {
735
736     Int_t i__1, i__2;
737     Float_t d__ [73] [73] ;
738     Int_t npad[201], i__, j, k;
739     Int_t id, id1, id2, iclb;
740
741 /*     using information from crude clustering, a refined algorithm for */
742 /*     clustering tried. At the moment using Gaussian fit with a number */
743 /*     of Gaussian functions. */
744         for (j = 1; j <= 72; ++j) {
745             for (k = 1; k <= 72; ++k) {
746              d__ [j] [k] = (*PPP).det2 [j] [k] ;
747             }
748         }
749
750     (*PPP).clno = 0;
751 /*     The No.of pads in each cluster = 200 */
752     for (i__ = 1; i__ <= 200; ++i__) {
753         npad[i__] = 0;
754     }
755     for (i__ = 1; i__ <= 5000; ++i__) {
756         (*PPP).ncl[i__] = 0;
757     }
758 /*     number of pads in each cluster determined. */
759     i__1 = (*PPP).incr;
760     for (i__ = 1; i__ <= i__1; ++i__) {
761         id = (*PPP).inford[1] [i__];
762         ++(*PPP).ncl[id];
763         Int_t nclid = (*PPP).ncl[id] ;
764     printf("id,ncl= %d %d\n",id,(*PPP).ncl[id]) ;
765     }
766     iclb = 0;
767 /*     number of clusters determined. */
768     for (i__ = 1; i__ <= 5000; ++i__) {
769         id = (*PPP).ncl[i__];
770         if (id != 0 ) {
771             ++iclb;
772             ++npad[id];
773         }
774     }
775     printf("iclb= %d \n",iclb) ;
776     id = 0;
777 /*     define coordinates of each pad in a cluster and call gausspara */
778 /*     to determine Gaussian parameters. ( should it be eliminated ? ) */
779     i__1 = iclb;
780     for (i__ = 1; i__ <= i__1; ++i__) {
781           i__2 = (*PPP).ncl[i__];
782           printf("(*PPP).ncl(i)= %d %d\n",i__,(*PPP).ncl[i__]) ;
783         for (j = 1; j <= i__2; ++j) {
784             ++id;
785             id1 = (*PPP).inford[2] [id];
786             id2 = (*PPP).inford[3] [id];
787             (*PPP).x[j] = (*PPP).coord[1] [id1] [id2];
788             (*PPP).y[j] = (*PPP).coord[2] [id1] [id2];
789             (*PPP).z__[j] = d__[id1] [id2];
790         }
791     printf(" I am in refclust calling ==gausspara_== \n") ;
792         gausspara_(i__2,PPP );
793     printf(" I am in refclust after ==gausspara_== \n") ;
794     }
795           printf("clno= %d \n",(*PPP).clno);
796     return ;
797 } /* refclust_ */
798
799 /* gausspara_ */ void gausspara_(Int_t n,PMD_dig_reco* PPP)
800 {
801     Int_t  i__1, i__2, i__3;
802     Float_t d__1, d__2;
803
804     Int_t   ihld, iord[501], idum, i__, j, k;
805     Int_t   itest;
806     Int_t   ig;
807     Float_t pi, xi, yi, xj, yj, rrr;
808     Float_t sum;
809
810 /*     single Gaussian of the form exp(-((x-x_0)^2+(y-y_0)^2)/r_0^2) */
811
812     printf(" I am in ==gausspara_== %d  \n" , n) ;
813     pi = 3.141593;
814     Int_t nn = n ;
815     if (nn == 1) {
816         ++(*PPP).clno;
817 /*     single pad cluster. Assumed fit to single Gaussian centered at the */
818 /*     center of the pad with Gaussian parameter r0 about 0.5 since */
819 /*     neighbouring pads don't fire. */
820         (*PPP).x0[1] = (*PPP).x[1];
821         (*PPP).y0[1] = (*PPP).y[1];
822         (*PPP).r0[1] = 0.5;
823 /* Computing 2nd power */
824         d__1 = (*PPP).r0[1];
825         (*PPP).h0[1] = (*PPP).z__[1] / pi / (d__1 * d__1);
826        printf("single pad ,clust no. = %d %d %f %f \n",n, (*PPP).clno, (*PPP).x0[1],(*PPP).y0[1] ) ;
827     } else if (nn == 2) {
828         ++((*PPP).clno);
829 /*     two pad cluster. Assumed fit to single Gaussian centered at the */
830 /*     weighted average of the centers of two pads with r_0 about 0.5. */
831 /*     Could be a pair of clusters but no way of differentiating the two. */
832 /*         print *,'two pads   ','clust no.  ', clno */
833        printf("two pad clusters  = %d \n",(*PPP).clno ) ;              
834 //       printf("z__(1),(2)  = %f %f \n",(*PPP).z__[1],(*PPP).z__[2] ) ; 
835        Int_t btest = ((*PPP).z__[1])+((*PPP).z__[2]) ;
836         if ( btest != 0 ) {
837         (*PPP).x0[1] = ((*PPP).z__[1] * (*PPP).x[1] + (*PPP).z__[2] * (*PPP).x[2]) / ((*PPP).z__[1] + (*PPP).z__[2]);
838         (*PPP).y0[1] = ((*PPP).z__[1] * (*PPP).y[1] + (*PPP).z__[2] * (*PPP).y[2]) / ((*PPP).z__[1] + (*PPP).z__[2]);
839         (*PPP).r0[1] = 0.5;
840         }
841 /* Computing 2nd power */
842         d__1 = (*PPP).r0[1] = 0.5;
843         (*PPP).h0[1] = ((*PPP).z__[1] + (*PPP).z__[2]) / pi / (d__1 * d__1);
844 //       printf("Two pad clusters  = %d %d %f %f \n",n, (*PPP).clno, (*PPP).x0[1],(*PPP).y0[1] ) ;
845     } else {
846 /*     guessing possible no. of Gaussians required. */
847 /*     assuming that if pads are separated by more than one unit, they */
848 /*     are possible candidates for Gaussian centers. */
849 /*     do this from larger values of energy ( z here ). */
850         i__1 = nn;
851         for (i__ = 1; i__ <= i__1; ++i__) {
852             iord[i__] = i__;
853         }
854         i__1 = nn;
855         for (i__ = 2; i__ <= i__1; ++i__) {
856             itest = 0;
857             ihld = iord[i__];
858             i__2 = i__ - 1;
859             for (j = 1; j <= i__2; ++j) {
860                 if (itest == 0 && (*PPP).z__[iord[j]] < (*PPP).z__[ihld]) {
861                     itest = 1;
862                     i__3 = j;
863                     for (k = i__ - 1; k >= i__3; --k) {
864                         iord[k+1] = iord[k];
865                     }
866                     iord[j] = ihld;
867                 }
868             }
869         }
870 /*     Gaussian centers determined. */
871         ig = 1;
872         (*PPP).icl[ig] = iord[ig];
873         i__1 = nn;
874         for (i__ = 2; i__ <= i__1; ++i__) {
875             xi = (*PPP).x[iord[i__]];
876             yi = (*PPP).y[iord[i__]];
877             idum = 0;
878             i__2 = ig;
879             for (j = 1; j <= i__2; ++j) {
880                 xj = (*PPP).x[(*PPP).icl[j]];
881                 yj = (*PPP).y[(*PPP).icl[j]];
882 /* Computing 2nd power */
883                 d__1 = xi - xj;
884 /* Computing 2nd power */
885                 d__2 = yi - yj;
886                 rrr = TMath::Sqrt(d__1 * d__1 + d__2 * d__2);
887                 if (1.01 < rrr) {
888                     ++idum;
889                 }
890             }
891             if (idum == ig) {
892                 ++ig;
893                 (*PPP).icl[ig] = iord[i__];
894             }
895         }
896         (*PPP).clno += ig;
897        printf("No. of pads = %d , No. of Gaussians = %f \n ",ig,(*PPP).clno ) ;
898         sum = 0.;
899         i__1 = ig;
900         for (i__ = 1; i__ <= i__1; ++i__) {
901             (*PPP).x0[i__] = (*PPP).x[(*PPP).icl[i__]];
902             (*PPP).y0[i__] = (*PPP).y[(*PPP).icl[i__]];
903             (*PPP).r0[i__] = .5;
904             sum += (*PPP).z__[(*PPP).icl[i__]];
905         }
906         i__1 = ig;
907         for (i__ = 1; i__ <= i__1; ++i__) {
908         if ( sum != 0 ) {
909             (*PPP).h0[i__] = (*PPP).z__[(*PPP).icl[i__]] / sum;
910         }
911         }
912     printf(" I am in ==gausspara_== Calling gausfit_ \n") ;
913         gausfit_(ig, nn, PPP );
914     printf(" I am in ==gausspara_ after gausfit_ \n") ;
915     }
916     return ;
917 } /* gausspara_ */
918
919 /*  gausfit_ */ void gausfit_(Int_t ng, Int_t nn, PMD_dig_reco* PPP )
920 {
921     Int_t  i__1;
922     Float_t d__1;
923
924     Int_t   iran;
925     Float_t aint;
926     Int_t   i__;
927     Float_t s;
928     Float_t  pi, ss, st;
929
930     pi = 3.141593;
931     printf(" ng=%d , nn=%d \n",ng,nn);
932     s = 0.;
933     i__1 = nn;
934     for (i__ = 1; i__ <= i__1; ++i__) {
935         s +=(*PPP).z__[i__];
936     }
937 /*      print *,s */
938      printf(" s=%d  \n",s);
939     i__1 = ng;
940     for (i__ = 1; i__ <= i__1; ++i__) {
941 /* Computing 2nd power */
942 //       printf("i__, X0,y0,r0 = %d %f %f %f \n",i__, (*PPP).x0[i__],(*PPP).y0[i__],(*PPP).r0[i__]) ;
943         d__1 = (*PPP).r0[i__];
944         (*PPP).h0[i__] = (*PPP).h0[i__] * s / pi / (d__1 * d__1);
945     }
946     st = 0.;
947     i__1 = nn;
948     for (i__ = 1; i__ <= i__1; ++i__) {
949         Int_t pbab = 1 ;
950         aintg_(ng, i__, pbab, aint, PPP);
951 /* Computing 2nd power */
952         d__1 = aint - (*PPP).z__[i__];
953         if( s > 0. ) st += d__1 * d__1 / s;
954     }
955 //   printf(" st=%d ",st);
956     for (iran = 1; iran <= 5000; ++iran) {
957         Float_t rand;
958         i__1 = ng;
959         for (i__ = 1; i__ <= i__1; ++i__) {
960             ranmar_(rand ,PPP) ;
961             (*PPP).xx[i__] = (*PPP).x0[i__] + (rand - .5) * 1.;
962             ranmar_(rand ,PPP) ;
963             (*PPP).yy[i__] = (*PPP).y0[i__] + (rand - .5) * 1.;
964             ranmar_(rand ,PPP) ;
965             (*PPP).rr[i__] = (*PPP).r0[i__] * (((rand -.5) * .2 )+ 1.);
966 /* Computing 2nd power */
967             d__1 = (*PPP).r0[i__] / (*PPP).rr[i__];
968             (*PPP).hh[i__] = (*PPP).h0[i__] * (d__1 * d__1);
969         }
970         ss = 0.;
971         i__1 = nn;
972         for (i__ = 1; i__ <= i__1; ++i__) {
973             Int_t pbab = 2 ;
974             aintg_(ng, i__, pbab, aint, PPP);
975 /* Computing 2nd power */
976             d__1 = aint - (*PPP).z__[i__];
977             if(s>0.0)ss += d__1 * d__1 / s;
978         }
979         if (ss < st) {
980             i__1 = ng;
981             for (i__ = 1; i__ <= i__1; ++i__) {
982                 (*PPP).x0[i__] = (*PPP).xx[i__];
983                 (*PPP).y0[i__] = (*PPP).yy[i__];
984                 (*PPP).r0[i__] = (*PPP).rr[i__];
985                 (*PPP).h0[i__] = (*PPP).hh[i__];
986             }
987             st = ss;
988     printf(" st=%d \n",st);
989         }
990     }
991 /*      print 1,((*PPP).x0(i),(*PPP).(*PPP).y0(i),(*PPP).r0(i),i=1,ng) */
992 /*      print *,st */
993             i__1 = ng;
994         for (i__ = 1; i__ <= i__1; ++i__) {
995 //             printf("x0(i),y0(i),r0(i) =%f %f %f \n",(*PPP).x0[i__1],(*PPP).y0[i__1],(*PPP).r0[i__1] ) ;
996         }
997     printf(" I am ending gausfit_ \n") ;
998     return ;
999 } /* gausfit_ */
1000
1001 /* aintg_ */void aintg_(Int_t ng, Int_t n, Int_t bab, Float_t aint, PMD_dig_reco* PPP )
1002 {
1003     Int_t i__1;
1004     Float_t  d__1, d__2, d__3;
1005
1006     Int_t i__, j;
1007     Float_t xp, yp, fun;
1008
1009     aint = 0.;
1010     for (i__ = 1; i__ <= 25; ++i__) {
1011         xp = (*PPP).pts[1] [i__] + (*PPP).x[n] ;
1012         yp = (*PPP).pts[2] [i__] + (*PPP).y[n] ;
1013 // printf("x[n],y[n],pts[1] [i__],[2] [i__], xp,yp =%f %f %f %f %f %f\n",(*PPP).x[n],(*PPP).y[n],(*PPP).pts[1] [i__],(*PPP).pts[2] [i__],xp,yp) ;
1014         fun = 0.;
1015         i__1 = ng;
1016         for (j = 1; j <= i__1; ++j) {
1017 //printf("J,D__1,D__2,D__3= %d %f %f %f %f\n",j,(*PPP).xx[j],xp,(*PPP).yy[j],(*PPP).rr[j]) ;
1018 /* Computing 2nd power */
1019             d__1 = xp - (*PPP).xx[j];
1020 /* Computing 2nd power */
1021             d__2 = yp - (*PPP).yy[j];
1022 /* Computing 2nd power */
1023             if( bab == 1) d__3 = (*PPP).r0[j];
1024             if( bab == 2) d__3 = (*PPP).rr[j];
1025 // printf("J,D__1,D__2,D__3= %d %f %f %f \n",j,d__1,d__2,d__3) ;
1026             fun += TMath::exp(-(d__1 * d__1 + d__2 * d__2) / (d__3 * d__3)) * (*PPP).hh[j];
1027         }
1028         aint += fun * (*PPP).wts[i__];
1029 //  printf("FUN,aint= %f %f \n",fun,aint) ;
1030     }
1031     return ;
1032 } /* aintg_ */
1033 /* ranmar() >*/void ranmar_(Float_t random , PMD_dig_reco* PPP)
1034 {
1035     /* Initialized data */
1036
1037     Int_t i__ = 97;
1038     Int_t j   = 33;
1039     Float_t    uni;
1040
1041 /* ===========================================C==========================C
1042  */
1043 /* Universal random number generator proposed by Marsaglia and Zaman */
1044 /* in report FSU-SCRI-87-50 */
1045     uni = (*PPP).u[i__ - 1] - (*PPP).u[j - 1];
1046     if (uni < 0.) {
1047         uni += 1.;
1048     }
1049     (*PPP).u[i__ ] = uni;
1050     --i__;
1051     if (i__ == 0) { i__ = 97; }
1052     --j;
1053     if (j   == 0) { j   = 97; }
1054     (*PPP).c__ -= (*PPP).cd;
1055     if ((*PPP).c__ < 0.) { (*PPP).c__ += (*PPP).cm; }
1056     uni -= (*PPP).c__;
1057     if (uni <= 0.) { uni += 1.; }
1058     random = uni;
1059 } /* ranmar_ */
1060