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