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