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