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