PMD/AliPMDDigitizer.cxx

   1 /***************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //-----------------------------------------------------//
  16 //                                                     //
  17 //  Source File : PMDDigitizer.cxx, Version 00         //
  18 //                                                     //
  19 //  Date   : September 20 2002                         //
  20 //                                                     //
  21 //-----------------------------------------------------//
  22
  23 #include <Riostream.h>
  24 #include <TBRIK.h>
  25 #include <TNode.h>
  26 #include <TTree.h>
  27 #include <TGeometry.h>
  28 #include <TObjArray.h>
  29 #include <TClonesArray.h>
  30 #include <TFile.h>
  31 #include <TNtuple.h>
  32 #include <TParticle.h>
  33
  34 #include "AliRun.h"
  35 #include "AliPMD.h"
  36 #include "AliHit.h"
  37 #include "AliDetector.h"
  38 #include "AliRunLoader.h"
  39 #include "AliLoader.h"
  40 #include "AliConfig.h"
  41 #include "AliMagF.h"
  42 #include "AliRunDigitizer.h"
  43 #include "AliHeader.h"
  44
  45 #include "AliPMDcell.h"
  46 #include "AliPMDsdigit.h"
  47 #include "AliPMDdigit.h"
  48 #include "AliPMDDigitizer.h"
  49 #include "AliPMDClustering.h"
  50 #include "AliMC.h"
  51
  52 ClassImp(AliPMDDigitizer)
  53 //
  54 // Constructor
  55 //
  56 AliPMDDigitizer::AliPMDDigitizer()
  57 {
  58   // Default Constructor
  59   //
  60   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
  61   fNsdigit = 0;
  62   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
  63   fNdigit = 0;
  64
  65   for (Int_t i = 0; i < fgkTotUM; i++)
  66     {
  67       for (Int_t j = 0; j < fgkRow; j++)
  68         {
  69           for (Int_t k = 0; k < fgkCol; k++)
  70             {
  71               fCPV[i][j][k] = 0.;
  72               fPRE[i][j][k] = 0.;
  73             }
  74         }
  75     }
  76
  77   if (!fCell) fCell = new TObjArray();
  78
  79   fZPos = 361.5; // in units of cm, This is the default position of PMD
  80 }
  81 AliPMDDigitizer::~AliPMDDigitizer()
  82 {
  83   // Default Destructor
  84   //
  85   delete fSDigits;
  86   delete fDigits;
  87   delete fCell;
  88 }
  89 //
  90 // Member functions
  91 //
  92 void AliPMDDigitizer::OpengAliceFile(Char_t *file, Option_t *option)
  93 {
  94   // Loads galice.root file and corresponding header, kinematics
  95   // hits and sdigits or digits depending on the option
  96   //
  97   fRunLoader = AliRunLoader::Open(file,AliConfig::fgkDefaultEventFolderName,
  98                                   "UPDATE");
  99
 100   if (!fRunLoader)
 101    {
 102      Error("Open","Can not open session for file %s.",file);
 103    }
 104
 105   fRunLoader->LoadgAlice();
 106   fRunLoader->LoadHeader();
 107   fRunLoader->LoadKinematics();
 108
 109   gAlice = fRunLoader->GetAliRun();
 110
 111   if (gAlice)
 112     {
 113       printf("<AliPMDdigitizer::Open> ");
 114       printf("AliRun object found on file.\n");
 115     }
 116   else
 117     {
 118       printf("<AliPMDdigitizer::Open> ");
 119       printf("Could not find AliRun object.\n");
 120     }
 121
 122   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
 123   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 124   if (fPMDLoader == 0x0)
 125     {
 126       cerr<<"Hits2Digits : Can not find PMD or PMDLoader\n";
 127     }
 128
 129   const char *cHS = strstr(option,"HS");
 130   const char *cHD = strstr(option,"HD");
 131   const char *cSD = strstr(option,"SD");
 132
 133   if (cHS)
 134     {
 135       fPMDLoader->LoadHits("READ");
 136       fPMDLoader->LoadSDigits("recreate");
 137     }
 138   else if (cHD)
 139     {
 140       fPMDLoader->LoadHits("READ");
 141       fPMDLoader->LoadDigits("recreate");
 142     }
 143   else if (cSD)
 144     {
 145       fPMDLoader->LoadSDigits("READ");
 146       fPMDLoader->LoadDigits("recreate");
 147     }
 148
 149 }
 150 void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
 151 {
 152   // This reads the PMD Hits tree and assigns the right track number
 153   // to a cell and stores in the summable digits tree
 154   //
 155   // cout << " -------- Beginning of Hits2SDigits ----------- " << endl;
 156
 157   const Int_t kPi0 = 111;
 158   const Int_t kGamma = 22;
 159   Int_t npmd;
 160   Int_t trackno;
 161   Int_t smnumber;
 162   Int_t trackpid;
 163   Int_t mtrackno;
 164   Int_t mtrackpid;
 165
 166   Float_t xPos, yPos, zPos;
 167   Int_t xpad = -1, ypad = -1;
 168   Float_t edep;
 169   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 170
 171
 172   ResetSDigit();
 173
 174   printf("Event Number =  %d \n",ievt);
 175   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
 176   printf("Number of Particles = %d \n", nparticles);
 177   fRunLoader->GetEvent(ievt);
 178   fPArray = gAlice->GetMCApp()->Particles();
 179   // ------------------------------------------------------- //
 180   // Pointer to specific detector hits.
 181   // Get pointers to Alice detectors and Hits containers
 182
 183   fTreeH = fPMDLoader->TreeH();
 184
 185   Int_t ntracks    = (Int_t) fTreeH->GetEntries();
 186   printf("Number of Tracks in the TreeH = %d \n", ntracks);
 187
 188   fTreeS = fPMDLoader->TreeS();
 189   if (fTreeS == 0x0)
 190     {
 191       fPMDLoader->MakeTree("S");
 192       fTreeS = fPMDLoader->TreeS();
 193     }
 194   Int_t bufsize = 16000;
 195   fTreeS->Branch("PMDSDigit", &fSDigits, bufsize);
 196
 197   if (fPMD) fHits   = fPMD->Hits();
 198
 199   // Start loop on tracks in the hits containers
 200
 201   for (Int_t track=0; track<ntracks;track++)
 202     {
 203       gAlice->ResetHits();
 204       fTreeH->GetEvent(track);
 205       if (fPMD)
 206         {
 207           npmd = fHits->GetEntriesFast();
 208           cout << " npmd = " << npmd << endl;
 209           for (int ipmd = 0; ipmd < npmd; ipmd++)
 210             {
 211               fPMDHit = (AliPMDhit*) fHits->UncheckedAt(ipmd);
 212               trackno = fPMDHit->GetTrack();
 213
 214               //  get kinematics of the particles
 215
 216               fParticle = gAlice->GetMCApp()->Particle(trackno);
 217               trackpid  = fParticle->GetPdgCode();
 218
 219               cout << " trackno = " << trackno << " track = " << track << endl;
 220               //cout << " trackpid = " << trackpid << endl;
 221               Int_t igatr = -999;
 222               Int_t ichtr = -999;
 223               Int_t igapid = -999;
 224               Int_t imo;
 225               Int_t igen = 0;
 226               Int_t idmo = -999;
 227
 228               TParticle*  mparticle = fParticle;
 229               Int_t jjnu = mparticle->GetFirstMother();
 230               cout << " Coming to this step - step 0 jjnu = " << jjnu << endl;
 231
 232               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
 233               if (mparticle->GetFirstMother() == -1)
 234                 {
 235                   tracknoOld  = trackno;
 236                   trackpidOld = trackpid;
 237                   statusOld   = -1;
 238                 }
 239               Int_t igstatus = 0;
 240               while((imo = mparticle->GetFirstMother()) >= 0)
 241                 {
 242                   igen++;
 243
 244                   mparticle =  gAlice->GetMCApp()->Particle(imo);
 245                   idmo = mparticle->GetPdgCode();
 246
 247                   vx = mparticle->Vx();
 248                   vy = mparticle->Vy();
 249                   vz = mparticle->Vz();
 250
 251                   //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 252                   //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 253                   if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
 254                     {
 255                       igatr = imo;
 256                       igapid = idmo;
 257                       igstatus = 1;
 258                     }
 259                   if(igstatus == 0)
 260                     {
 261                       if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 262                         {
 263                           igatr = imo;
 264                           igapid = idmo;
 265                         }
 266                     }
 267                   ichtr = imo;
 268                 }
 269
 270               if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 271                 {
 272                   mtrackno = igatr;
 273                   mtrackpid = igapid;
 274                 }
 275               else
 276                 {
 277                   mtrackno  = ichtr;
 278                   mtrackpid = idmo;
 279                 }
 280               if (statusOld == -1)
 281                 {
 282                   mtrackno  = tracknoOld;
 283                   mtrackpid = trackpidOld;
 284                 }
 285               xPos = fPMDHit->X();
 286               yPos = fPMDHit->Y();
 287               zPos = fPMDHit->Z();
 288               edep       = fPMDHit->fEnergy;
 289               Int_t vol1 = fPMDHit->fVolume[1]; // Column
 290               Int_t vol2 = fPMDHit->fVolume[2]; // Row
 291               Int_t vol3 = fPMDHit->fVolume[3]; // UnitModule
 292               Int_t vol6 = fPMDHit->fVolume[6]; // SuperModule
 293               // -----------------------------------------//
 294               // For Super Module 1 & 2                   //
 295               //  nrow = 96, ncol = 48                    //
 296               // For Super Module 3 & 4                   //
 297               //  nrow = 48, ncol = 96                    //
 298               // -----------------------------------------//
 299
 300               smnumber = (vol6-1)*6 + vol3;
 301
 302               if (vol6 == 1 || vol6 == 2)
 303                 {
 304                   xpad = vol1;
 305                   ypad = vol2;
 306                 }
 307               else if (vol6 == 3 || vol6 == 4)
 308                 {
 309                   xpad = vol2;
 310                   ypad = vol1;
 311                 }
 312
 313               //cout << "zpos = " << zPos << " edep = " << edep << endl;
 314
 315               Float_t zposition = TMath::Abs(zPos);
 316               if (zposition < fZPos)
 317                 {
 318                   // CPV
 319                   fDetNo = 1;
 320                 }
 321               else if (zposition > fZPos)
 322                 {
 323                   // PMD
 324                   fDetNo = 0;
 325                 }
 326               Int_t smn = smnumber - 1;
 327               Int_t ixx = xpad     - 1;
 328               Int_t iyy = ypad     - 1;
 329               if (fDetNo == 0)
 330                 {
 331                   fPRE[smn][ixx][iyy] += edep;
 332                   fPRECounter[smn][ixx][iyy]++;
 333
 334                   fPMDcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 335
 336                   fCell->Add(fPMDcell);
 337                 }
 338               else if(fDetNo == 1)
 339                 {
 340                   fCPV[smn][ixx][iyy] += edep;
 341                   fCPVTrackNo[smn][ixx][iyy] = mtrackno;
 342                 }
 343             }
 344         }
 345     } // Track Loop ended
 346
 347   TrackAssignment2Cell();
 348   ResetCell();
 349
 350   Float_t deltaE      = 0.;
 351   Int_t   detno       = 0;
 352   Int_t   trno        = -1;
 353   Int_t   cellno      = 0;
 354
 355   for (Int_t idet = 0; idet < 2; idet++)
 356     {
 357       for (Int_t ism = 0; ism < fgkTotUM; ism++)
 358         {
 359           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
 360             {
 361               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
 362                 {
 363                   cellno = jrow*fgkCol + kcol;
 364                   if (idet == 0)
 365                     {
 366                       deltaE = fPRE[ism][jrow][kcol];
 367                       trno   = fPRETrackNo[ism][jrow][kcol];
 368                       detno = 0;
 369                     }
 370                   else if (idet == 1)
 371                     {
 372                       deltaE = fCPV[ism][jrow][kcol];
 373                       trno   = fCPVTrackNo[ism][jrow][kcol];
 374                       detno = 1;
 375                     }
 376                   if (deltaE > 0.)
 377                     {
 378                       AddSDigit(trno,detno,ism,cellno,deltaE);
 379                     }
 380                 }
 381             }
 382           fTreeS->Fill();
 383           ResetSDigit();
 384         }
 385     }
 386   fPMDLoader->WriteSDigits("OVERWRITE");
 387   ResetCellADC();
 388
 389   //  cout << " -------- End of Hits2SDigit ----------- " << endl;
 390 }
 391
 392 void AliPMDDigitizer::Hits2Digits(Int_t ievt)
 393 {
 394   // This reads the PMD Hits tree and assigns the right track number
 395   // to a cell and stores in the digits tree
 396   //
 397   const Int_t kPi0 = 111;
 398   const Int_t kGamma = 22;
 399   Int_t npmd;
 400   Int_t trackno;
 401   Int_t smnumber;
 402   Int_t trackpid;
 403   Int_t mtrackno;
 404   Int_t mtrackpid;
 405
 406   Float_t xPos, yPos, zPos;
 407   Int_t xpad = -1, ypad = -1;
 408   Float_t edep;
 409   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 410
 411
 412   ResetDigit();
 413
 414   printf("Event Number =  %d \n",ievt);
 415
 416   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
 417   printf("Number of Particles = %d \n", nparticles);
 418   fRunLoader->GetEvent(ievt);
 419   fPArray = gAlice->GetMCApp()->Particles();
 420   // ------------------------------------------------------- //
 421   // Pointer to specific detector hits.
 422   // Get pointers to Alice detectors and Hits containers
 423
 424   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
 425   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 426
 427   if (fPMDLoader == 0x0)
 428     {
 429       cerr<<"Hits2Digits method : Can not find PMD or PMDLoader\n";
 430     }
 431   fTreeH = fPMDLoader->TreeH();
 432   Int_t ntracks    = (Int_t) fTreeH->GetEntries();
 433   printf("Number of Tracks in the TreeH = %d \n", ntracks);
 434   fPMDLoader->LoadDigits("recreate");
 435   fTreeD = fPMDLoader->TreeD();
 436   if (fTreeD == 0x0)
 437     {
 438       fPMDLoader->MakeTree("D");
 439       fTreeD = fPMDLoader->TreeD();
 440     }
 441   Int_t bufsize = 16000;
 442   fTreeD->Branch("PMDDigit", &fDigits, bufsize);
 443
 444   if (fPMD) fHits   = fPMD->Hits();
 445
 446   // Start loop on tracks in the hits containers
 447
 448   for (Int_t track=0; track<ntracks;track++)
 449     {
 450       gAlice->ResetHits();
 451       fTreeH->GetEvent(track);
 452
 453       if (fPMD)
 454         {
 455           npmd = fHits->GetEntriesFast();
 456           for (int ipmd = 0; ipmd < npmd; ipmd++)
 457             {
 458               fPMDHit = (AliPMDhit*) fHits->UncheckedAt(ipmd);
 459               trackno = fPMDHit->GetTrack();
 460
 461               //  get kinematics of the particles
 462
 463               fParticle = gAlice->GetMCApp()->Particle(trackno);
 464               trackpid  = fParticle->GetPdgCode();
 465
 466               Int_t igatr = -999;
 467               Int_t ichtr = -999;
 468               Int_t igapid = -999;
 469               Int_t imo;
 470               Int_t igen = 0;
 471               Int_t idmo = -999;
 472
 473               TParticle*  mparticle = fParticle;
 474               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
 475               if (mparticle->GetFirstMother() == -1)
 476                 {
 477                   tracknoOld  = trackno;
 478                   trackpidOld = trackpid;
 479                   statusOld   = -1;
 480                 }
 481
 482               Int_t igstatus = 0;
 483               while((imo = mparticle->GetFirstMother()) >= 0)
 484                 {
 485                   igen++;
 486
 487                   mparticle =  gAlice->GetMCApp()->Particle(imo);
 488                   idmo = mparticle->GetPdgCode();
 489
 490                   vx = mparticle->Vx();
 491                   vy = mparticle->Vy();
 492                   vz = mparticle->Vz();
 493
 494                   //printf("==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 495                   //fprintf(fpw1,"==> Mother ID %5d %5d %5d Vertex: %13.3f %13.3f %13.3f\n", igen, imo, idmo, vx, vy, vz);
 496                   if ((idmo == kGamma || idmo == -11 || idmo == 11) && vx == 0. && vy == 0. && vz == 0.)
 497                     {
 498                       igatr = imo;
 499                       igapid = idmo;
 500                       igstatus = 1;
 501                     }
 502                   if(igstatus == 0)
 503                     {
 504                       if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 505                         {
 506                           igatr = imo;
 507                           igapid = idmo;
 508                         }
 509                     }
 510                   ichtr = imo;
 511                 }
 512
 513               if (idmo == kPi0 && vx == 0. && vy == 0. && vz == 0.)
 514                 {
 515                   mtrackno = igatr;
 516                   mtrackpid = igapid;
 517                 }
 518               else
 519                 {
 520                   mtrackno  = ichtr;
 521                   mtrackpid = idmo;
 522                 }
 523               if (statusOld == -1)
 524                 {
 525                   mtrackno  = tracknoOld;
 526                   mtrackpid = trackpidOld;
 527                 }
 528
 529               xPos = fPMDHit->X();
 530               yPos = fPMDHit->Y();
 531               zPos = fPMDHit->Z();
 532               edep       = fPMDHit->fEnergy;
 533               Int_t vol1 = fPMDHit->fVolume[1]; // Column
 534               Int_t vol2 = fPMDHit->fVolume[2]; // Row
 535               Int_t vol3 = fPMDHit->fVolume[3]; // UnitModule
 536               Int_t vol6 = fPMDHit->fVolume[6]; // SuperModule
 537
 538               // -----------------------------------------//
 539               // For Super Module 1 & 2                   //
 540               //  nrow = 96, ncol = 48                    //
 541               // For Super Module 3 & 4                   //
 542               //  nrow = 48, ncol = 96                    //
 543               // -----------------------------------------//
 544
 545               smnumber = (vol6-1)*6 + vol3;
 546
 547               if (vol6 == 1 || vol6 == 2)
 548                 {
 549                   xpad = vol1;
 550                   ypad = vol2;
 551                 }
 552               else if (vol6 == 3 || vol6 == 4)
 553                 {
 554                   xpad = vol2;
 555                   ypad = vol1;
 556                 }
 557
 558               //cout << "-zpos = " << -zPos << endl;
 559
 560               Float_t zposition = TMath::Abs(zPos);
 561
 562               if (zposition < fZPos)
 563                 {
 564                   // CPV
 565                   fDetNo = 1;
 566                 }
 567               else if (zposition > fZPos)
 568                 {
 569                   // PMD
 570                   fDetNo = 0;
 571                 }
 572
 573               Int_t smn = smnumber - 1;
 574               Int_t ixx = xpad     - 1;
 575               Int_t iyy = ypad     - 1;
 576               if (fDetNo == 0)
 577                 {
 578                   fPRE[smn][ixx][iyy] += edep;
 579                   fPRECounter[smn][ixx][iyy]++;
 580
 581                   fPMDcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 582
 583                   fCell->Add(fPMDcell);
 584                 }
 585               else if(fDetNo == 1)
 586                 {
 587                   fCPV[smn][ixx][iyy] += edep;
 588                   fCPVTrackNo[smn][ixx][iyy] = mtrackno;
 589                 }
 590             }
 591         }
 592     } // Track Loop ended
 593
 594   TrackAssignment2Cell();
 595   ResetCell();
 596
 597   Float_t deltaE = 0.;
 598   Int_t detno = 0;
 599   Int_t trno = 1;
 600   Int_t cellno;
 601
 602   for (Int_t idet = 0; idet < 2; idet++)
 603     {
 604       for (Int_t ism = 0; ism < fgkTotUM; ism++)
 605         {
 606           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
 607             {
 608               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
 609                 {
 610                   cellno = jrow*fgkCol + kcol;
 611                   if (idet == 0)
 612                     {
 613                       deltaE = fPRE[ism][jrow][kcol];
 614                       trno   = fPRETrackNo[ism][jrow][kcol];
 615                       detno = 0;
 616                     }
 617                   else if (idet == 1)
 618                     {
 619                       deltaE = fCPV[ism][jrow][kcol];
 620                       trno   = fCPVTrackNo[ism][jrow][kcol];
 621                       detno = 1;
 622                     }
 623                   if (deltaE > 0.)
 624                     {
 625                       AddDigit(trno,detno,ism,cellno,deltaE);
 626                     }
 627                 } // column loop
 628             } // row    loop
 629         } // supermodule loop
 630       fTreeD->Fill();
 631       ResetDigit();
 632     } // detector loop
 633
 634   fPMDLoader->WriteDigits("OVERWRITE");
 635   ResetCellADC();
 636
 637   //  cout << " -------- End of Hits2Digit ----------- " << endl;
 638 }
 639
 640
 641 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
 642 {
 643   // This reads the PMD sdigits tree and converts energy deposition
 644   // in a cell to ADC and stores in the digits tree
 645   //
 646   //  cout << " -------- Beginning of SDigits2Digit ----------- " << endl;
 647   fRunLoader->GetEvent(ievt);
 648
 649   fTreeS = fPMDLoader->TreeS();
 650   AliPMDsdigit  *pmdsdigit;
 651   TBranch *branch = fTreeS->GetBranch("PMDSDigit");
 652   branch->SetAddress(&fSDigits);
 653
 654   fTreeD = fPMDLoader->TreeD();
 655   if (fTreeD == 0x0)
 656     {
 657       fPMDLoader->MakeTree("D");
 658       fTreeD = fPMDLoader->TreeD();
 659     }
 660   Int_t bufsize = 16000;
 661   fTreeD->Branch("PMDDigit", &fDigits, bufsize);
 662
 663   Int_t   trno, det, smn;
 664   Int_t   cellno;
 665   Float_t edep, adc;
 666
 667   Int_t nmodules = (Int_t) fTreeS->GetEntries();
 668
 669   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 670     {
 671       fTreeS->GetEntry(imodule);
 672       Int_t nentries = fSDigits->GetLast();
 673       //cout << " nentries = " << nentries << endl;
 674       for (Int_t ient = 0; ient < nentries+1; ient++)
 675         {
 676           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
 677           trno   = pmdsdigit->GetTrackNumber();
 678           det    = pmdsdigit->GetDetector();
 679           smn    = pmdsdigit->GetSMNumber();
 680           cellno = pmdsdigit->GetCellNumber();
 681           edep   = pmdsdigit->GetCellEdep();
 682
 683           MeV2ADC(edep,adc);
 684           AddDigit(trno,det,smn,cellno,adc);
 685         }
 686       fTreeD->Fill();
 687       ResetDigit();
 688     }
 689   fPMDLoader->WriteDigits("OVERWRITE");
 690   //  cout << " -------- End of SDigits2Digit ----------- " << endl;
 691 }
 692
 693 void AliPMDDigitizer::TrackAssignment2Cell()
 694 {
 695   //
 696   // This block assigns the cell id when there are
 697   // multiple tracks in a cell according to the
 698   // energy deposition
 699   //
 700   Bool_t jsort = false;
 701
 702   Int_t i, j, k;
 703
 704   Float_t *fracEdp;
 705   Float_t *trEdp;
 706   Int_t *status1;
 707   Int_t *status2;
 708   Int_t *trnarray;
 709   Int_t   ****pmdTrack;
 710   Float_t ****pmdEdep;
 711
 712   pmdTrack = new Int_t ***[fgkTotUM];
 713   pmdEdep  = new Float_t ***[fgkTotUM];
 714   for (i=0; i<fgkTotUM; i++)
 715     {
 716       pmdTrack[i] = new Int_t **[fgkRow];
 717       pmdEdep[i]  = new Float_t **[fgkRow];
 718     }
 719
 720   for (i = 0; i < fgkTotUM; i++)
 721     {
 722       for (j = 0; j < fgkRow; j++)
 723         {
 724           pmdTrack[i][j] = new Int_t *[fgkCol];
 725           pmdEdep[i][j]  = new Float_t *[fgkCol];
 726         }
 727     }
 728
 729   for (i = 0; i < fgkTotUM; i++)
 730     {
 731       for (j = 0; j < fgkRow; j++)
 732         {
 733           for (k = 0; k < fgkCol; k++)
 734             {
 735               Int_t nn = fPRECounter[i][j][k];
 736               if(nn > 0)
 737                 {
 738                   pmdTrack[i][j][k] = new Int_t[nn];
 739                   pmdEdep[i][j][k] = new Float_t[nn];
 740                 }
 741               else
 742                 {
 743                   nn = 1;
 744                   pmdTrack[i][j][k] = new Int_t[nn];
 745                   pmdEdep[i][j][k] = new Float_t[nn];
 746                 }
 747               fPRECounter[i][j][k] = 0;
 748             }
 749         }
 750     }
 751
 752
 753   Int_t nentries = fCell->GetEntries();
 754
 755   Int_t   mtrackno, ism, ixp, iyp;
 756   Float_t edep;
 757
 758   for (i = 0; i < nentries; i++)
 759     {
 760       fPMDcell = (AliPMDcell*)fCell->UncheckedAt(i);
 761
 762       mtrackno = fPMDcell->GetTrackNumber();
 763       ism      = fPMDcell->GetSMNumber();
 764       ixp      = fPMDcell->GetX();
 765       iyp      = fPMDcell->GetY();
 766       edep     = fPMDcell->GetEdep();
 767       Int_t nn = fPRECounter[ism][ixp][iyp];
 768       //      cout << " nn = " << nn << endl;
 769       pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
 770       pmdEdep[ism][ixp][iyp][nn] = edep;
 771       fPRECounter[ism][ixp][iyp]++;
 772     }
 773
 774   Int_t iz, il;
 775   Int_t im, ix, iy;
 776   Int_t nn;
 777
 778   for (im=0; im<fgkTotUM; im++)
 779     {
 780       for (ix=0; ix<fgkRow; ix++)
 781         {
 782           for (iy=0; iy<fgkCol; iy++)
 783             {
 784               nn = fPRECounter[im][ix][iy];
 785               if (nn > 1)
 786                 {
 787                   // This block handles if a cell is fired
 788                   // many times by many tracks
 789                   status1  = new Int_t[nn];
 790                   status2  = new Int_t[nn];
 791                   trnarray = new Int_t[nn];
 792                   for (iz = 0; iz < nn; iz++)
 793                     {
 794                       status1[iz] = pmdTrack[im][ix][iy][iz];
 795                     }
 796                   TMath::Sort(nn,status1,status2,jsort);
 797                   Int_t trackOld = -99999;
 798                   Int_t track, trCount = 0;
 799                   for (iz = 0; iz < nn; iz++)
 800                     {
 801                       track = status1[status2[iz]];
 802                       if (trackOld != track)
 803                         {
 804                           trnarray[trCount] = track;
 805                           trCount++;
 806                         }
 807                       trackOld = track;
 808                     }
 809                   delete status1;
 810                   delete status2;
 811                   Float_t totEdp = 0.;
 812                   trEdp = new Float_t[trCount];
 813                   fracEdp = new Float_t[trCount];
 814                   for (il = 0; il < trCount; il++)
 815                     {
 816                       trEdp[il] = 0.;
 817                       track = trnarray[il];
 818                       for (iz = 0; iz < nn; iz++)
 819                         {
 820                           if (track == pmdTrack[im][ix][iy][iz])
 821                             {
 822                               trEdp[il] += pmdEdep[im][ix][iy][iz];
 823                             }
 824                         }
 825                       totEdp += trEdp[il];
 826                     }
 827                   Int_t ilOld = 0;
 828                   Float_t fracOld = 0.;
 829
 830                   for (il = 0; il < trCount; il++)
 831                     {
 832                       fracEdp[il] = trEdp[il]/totEdp;
 833                       if (fracOld < fracEdp[il])
 834                         {
 835                           fracOld = fracEdp[il];
 836                           ilOld = il;
 837                         }
 838                     }
 839                   fPRETrackNo[im][ix][iy] = trnarray[ilOld];
 840                   delete fracEdp;
 841                   delete trEdp;
 842                   delete trnarray;
 843                 }
 844               else if (nn == 1)
 845                 {
 846                   // This only handles if a cell is fired
 847                   // by only one track
 848
 849                   fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
 850
 851                 }
 852               else if (nn ==0)
 853                 {
 854                   // This is if no cell is fired
 855                   fPRETrackNo[im][ix][iy] = -999;
 856                 }
 857             } // end of iy
 858         } // end of ix
 859     } // end of im
 860
 861   // Delete all the pointers
 862
 863   for (i = 0; i < fgkTotUM; i++)
 864     {
 865       for (j = 0; j < fgkRow; j++)
 866         {
 867           for (k = 0; k < fgkCol; k++)
 868             {
 869               delete [] pmdTrack[i][j][k];
 870               delete [] pmdEdep[i][j][k];
 871             }
 872         }
 873     }
 874
 875   for (i = 0; i < fgkTotUM; i++)
 876     {
 877       for (j = 0; j < fgkRow; j++)
 878         {
 879           delete [] pmdTrack[i][j];
 880           delete [] pmdEdep[i][j];
 881         }
 882     }
 883
 884   for (i = 0; i < fgkTotUM; i++)
 885     {
 886       delete [] pmdTrack[i];
 887       delete [] pmdEdep[i];
 888     }
 889   delete pmdTrack;
 890   delete pmdEdep;
 891   //
 892   // End of the cell id assignment
 893   //
 894 }
 895
 896
 897 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc)
 898 {
 899   // This converts the simulated edep to ADC according to the
 900   // Test Beam Data
 901   // To be done
 902   //
 903   adc = mev*1.;
 904 }
 905 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t det, Int_t smnumber,
 906   Int_t cellnumber, Float_t adc)
 907 {
 908   // Add SDigit
 909   //
 910   TClonesArray &lsdigits = *fSDigits;
 911   AliPMDsdigit *newcell;
 912   newcell = new AliPMDsdigit(trnumber,det,smnumber,cellnumber,adc);
 913   new(lsdigits[fNsdigit++]) AliPMDsdigit(newcell);
 914   delete newcell;
 915 }
 916
 917 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t det, Int_t smnumber,
 918   Int_t cellnumber, Float_t adc)
 919 {
 920   // Add Digit
 921   //
 922   TClonesArray &ldigits = *fDigits;
 923   AliPMDdigit *newcell;
 924   newcell = new AliPMDdigit(trnumber,det,smnumber,cellnumber,adc);
 925   new(ldigits[fNdigit++]) AliPMDdigit(newcell);
 926   delete newcell;
 927 }
 928
 929 void AliPMDDigitizer::SetZPosition(Float_t zpos)
 930 {
 931   fZPos = zpos;
 932 }
 933 Float_t AliPMDDigitizer::GetZPosition() const
 934 {
 935   return fZPos;
 936 }
 937
 938 void AliPMDDigitizer::ResetCell()
 939 {
 940   // clears the cell array and also the counter
 941   //  for each cell
 942   //
 943   fCell->Clear();
 944   for (Int_t i = 0; i < fgkTotUM; i++)
 945     {
 946       for (Int_t j = 0; j < fgkRow; j++)
 947         {
 948           for (Int_t k = 0; k < fgkCol; k++)
 949             {
 950               fPRECounter[i][j][k] = 0;
 951             }
 952         }
 953     }
 954 }
 955 void AliPMDDigitizer::ResetSDigit()
 956 {
 957   // Clears SDigits
 958   fNsdigit = 0;
 959   if (fSDigits) fSDigits->Clear();
 960 }
 961 void AliPMDDigitizer::ResetDigit()
 962 {
 963   // Clears Digits
 964   fNdigit = 0;
 965   if (fDigits) fDigits->Clear();
 966 }
 967
 968 void AliPMDDigitizer::ResetCellADC()
 969 {
 970   // Clears individual cells edep
 971   for (Int_t i = 0; i < fgkTotUM; i++)
 972     {
 973       for (Int_t j = 0; j < fgkRow; j++)
 974         {
 975           for (Int_t k = 0; k < fgkCol; k++)
 976             {
 977               fCPV[i][j][k] = 0.;
 978               fPRE[i][j][k] = 0.;
 979             }
 980         }
 981     }
 982 }
 983
 984 void AliPMDDigitizer::UnLoad(Option_t *option)
 985 {
 986   // Unloads all the root files
 987   //
 988   const char *cS = strstr(option,"S");
 989   const char *cD = strstr(option,"D");
 990
 991   fRunLoader->UnloadgAlice();
 992   fRunLoader->UnloadHeader();
 993   fRunLoader->UnloadKinematics();
 994
 995   if (cS)
 996     {
 997       fPMDLoader->UnloadHits();
 998     }
 999   if (cD)
1000     {
1001       fPMDLoader->UnloadHits();
1002       fPMDLoader->UnloadSDigits();
1003     }
1004 }