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 #include <TRandom.h>
  34
  35 #include "AliLog.h"
  36 #include "AliRun.h"
  37 #include "AliHit.h"
  38 #include "AliDetector.h"
  39 #include "AliRunLoader.h"
  40 #include "AliLoader.h"
  41 #include "AliConfig.h"
  42 #include "AliMagF.h"
  43 #include "AliRunDigitizer.h"
  44 #include "AliDigitizer.h"
  45 #include "AliHeader.h"
  46 #include "AliCDBManager.h"
  47 #include "AliCDBStorage.h"
  48 #include "AliCDBEntry.h"
  49 #include "AliMC.h"
  50
  51 #include "AliPMD.h"
  52 #include "AliPMDhit.h"
  53 #include "AliPMDcell.h"
  54 #include "AliPMDsdigit.h"
  55 #include "AliPMDdigit.h"
  56 #include "AliPMDCalibData.h"
  57 #include "AliPMDPedestal.h"
  58 #include "AliPMDDigitizer.h"
  59
  60
  61 ClassImp(AliPMDDigitizer)
  62
  63 AliPMDDigitizer::AliPMDDigitizer() :
  64   fRunLoader(0),
  65   fPMDHit(0),
  66   fPMD(0),
  67   fPMDLoader(0),
  68   fCalibGain(GetCalibGain()),
  69   fCalibPed(GetCalibPed()),
  70   fSDigits(0),
  71   fDigits(0),
  72   fCPVCell(0),
  73   fCell(0),
  74   fNsdigit(0),
  75   fNdigit(0),
  76   fDetNo(0),
  77   fZPos(361.5)   // in units of cm, default position of PMD
  78 {
  79   // Default Constructor
  80   //
  81   for (Int_t i = 0; i < fgkTotUM; i++)
  82     {
  83       for (Int_t j = 0; j < fgkRow; j++)
  84         {
  85           for (Int_t k = 0; k < fgkCol; k++)
  86             {
  87               fCPV[i][j][k]         = 0.;
  88               fPRE[i][j][k]         = 0.;
  89               fCPVCounter[i][j][k]  =  0;
  90               fPRECounter[i][j][k]  =  0;
  91               fCPVTrackNo[i][j][k]  = -1;
  92               fPRETrackNo[i][j][k]  = -1;
  93               fCPVTrackPid[i][j][k] = -1;
  94               fPRETrackPid[i][j][k] = -1;
  95             }
  96         }
  97     }
  98
  99 }
 100 //____________________________________________________________________________
 101 AliPMDDigitizer::AliPMDDigitizer(const AliPMDDigitizer& digitizer):
 102   AliDigitizer(digitizer),
 103   fRunLoader(0),
 104   fPMDHit(0),
 105   fPMD(0),
 106   fPMDLoader(0),
 107   fCalibGain(GetCalibGain()),
 108   fCalibPed(GetCalibPed()),
 109   fSDigits(0),
 110   fDigits(0),
 111   fCPVCell(0),
 112   fCell(0),
 113   fNsdigit(0),
 114   fNdigit(0),
 115   fDetNo(0),
 116   fZPos(361.5)   // in units of cm, default position of PMD
 117 {
 118   // copy constructor
 119   AliError("Copy constructor not allowed ");
 120
 121 }
 122 //____________________________________________________________________________
 123 AliPMDDigitizer & AliPMDDigitizer::operator=(const AliPMDDigitizer& /*digitizer*/)
 124 {
 125   // Assignment operator
 126   AliError("Assignement operator not allowed ");
 127
 128   return *this;
 129 }
 130 //____________________________________________________________________________
 131 AliPMDDigitizer::AliPMDDigitizer(AliRunDigitizer* manager):
 132   AliDigitizer(manager),
 133   fRunLoader(0),
 134   fPMDHit(0),
 135   fPMD(0),
 136   fPMDLoader(0),
 137   fCalibGain(GetCalibGain()),
 138   fCalibPed(GetCalibPed()),
 139   fSDigits(new TClonesArray("AliPMDsdigit", 1000)),
 140   fDigits(new TClonesArray("AliPMDdigit", 1000)),
 141   fCPVCell(0),
 142   fCell(0),
 143   fNsdigit(0),
 144   fNdigit(0),
 145   fDetNo(0),
 146   fZPos(361.5)// in units of cm, This is the default position of PMD
 147 {
 148   // ctor which should be used
 149
 150
 151   for (Int_t i = 0; i < fgkTotUM; i++)
 152     {
 153       for (Int_t j = 0; j < fgkRow; j++)
 154         {
 155           for (Int_t k = 0; k < fgkCol; k++)
 156             {
 157               fCPV[i][j][k]         = 0.;
 158               fPRE[i][j][k]         = 0.;
 159               fCPVCounter[i][j][k]  =  0;
 160               fPRECounter[i][j][k]  =  0;
 161               fCPVTrackNo[i][j][k]  = -1;
 162               fPRETrackNo[i][j][k]  = -1;
 163               fCPVTrackPid[i][j][k] = -1;
 164               fPRETrackPid[i][j][k] = -1;
 165             }
 166         }
 167     }
 168 }
 169
 170 //____________________________________________________________________________
 171 AliPMDDigitizer::~AliPMDDigitizer()
 172 {
 173   // Default Destructor
 174   //
 175   if (fSDigits) {
 176     fSDigits->Delete();
 177     delete fSDigits;
 178     fSDigits=0;
 179   }
 180   if (fDigits) {
 181     fDigits->Delete();
 182     delete fDigits;
 183     fDigits=0;
 184   }
 185   fCPVCell.Delete();
 186   fCell.Delete();
 187 }
 188 //
 189 // Member functions
 190 //
 191 //____________________________________________________________________________
 192 void AliPMDDigitizer::OpengAliceFile(const char *file, Option_t *option)
 193 {
 194   // Loads galice.root file and corresponding header, kinematics
 195   // hits and sdigits or digits depending on the option
 196   //
 197
 198   TString evfoldname = AliConfig::GetDefaultEventFolderName();
 199   fRunLoader = AliRunLoader::GetRunLoader(evfoldname);
 200   if (!fRunLoader)
 201       fRunLoader = AliRunLoader::Open(file,AliConfig::GetDefaultEventFolderName(), "UPDATE");
 202
 203   if (!fRunLoader)
 204    {
 205      AliError(Form("Can not open session for file %s.",file));
 206    }
 207
 208   const char *cHS = strstr(option,"HS");
 209   const char *cHD = strstr(option,"HD");
 210   const char *cSD = strstr(option,"SD");
 211
 212   if(cHS || cHD)
 213     {
 214       if (!fRunLoader->GetAliRun()) fRunLoader->LoadgAlice();
 215       if (!fRunLoader->TreeE()) fRunLoader->LoadHeader();
 216       if (!fRunLoader->TreeK()) fRunLoader->LoadKinematics();
 217
 218       gAlice = fRunLoader->GetAliRun();
 219
 220       if (gAlice)
 221         {
 222           AliDebug(1,"Alirun object found");
 223         }
 224       else
 225         {
 226           AliError("Could not found Alirun object");
 227         }
 228
 229       fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
 230     }
 231
 232   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 233   if (fPMDLoader == 0x0)
 234     {
 235       AliError("Can not find PMDLoader");
 236     }
 237
 238
 239   if (cHS)
 240     {
 241       fPMDLoader->LoadHits("READ");
 242       fPMDLoader->LoadSDigits("recreate");
 243     }
 244   else if (cHD)
 245     {
 246       fPMDLoader->LoadHits("READ");
 247       fPMDLoader->LoadDigits("recreate");
 248     }
 249   else if (cSD)
 250     {
 251       fPMDLoader->LoadSDigits("READ");
 252       fPMDLoader->LoadDigits("recreate");
 253     }
 254 }
 255 //____________________________________________________________________________
 256 void AliPMDDigitizer::Hits2SDigits(Int_t ievt)
 257 {
 258   // This reads the PMD Hits tree and assigns the right track number
 259   // to a cell and stores in the summable digits tree
 260   //
 261
 262   const Int_t kPi0 = 111;
 263   const Int_t kGamma = 22;
 264   Int_t npmd;
 265   Int_t trackno;
 266   Int_t smnumber;
 267   Int_t trackpid;
 268   Int_t mtrackno;
 269   Int_t mtrackpid;
 270
 271   Float_t xPos, yPos, zPos;
 272   Int_t xpad = -1, ypad = -1;
 273   Float_t edep;
 274   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 275
 276
 277   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
 278   ResetSDigit();
 279
 280   AliDebug(1,Form("Event Number = %d",ievt));
 281   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
 282   AliDebug(1,Form("Number of Particles = %d",nparticles));
 283
 284   //
 285
 286   fRunLoader->GetEvent(ievt);
 287   // ------------------------------------------------------- //
 288   // Pointer to specific detector hits.
 289   // Get pointers to Alice detectors and Hits containers
 290
 291   TTree* treeH = fPMDLoader->TreeH();
 292
 293   Int_t ntracks    = (Int_t) treeH->GetEntries();
 294   AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
 295   TTree* treeS = fPMDLoader->TreeS();
 296   if (treeS == 0x0)
 297     {
 298       fPMDLoader->MakeTree("S");
 299       treeS = fPMDLoader->TreeS();
 300     }
 301   Int_t bufsize = 16000;
 302   treeS->Branch("PMDSDigit", &fSDigits, bufsize);
 303
 304   TClonesArray* hits = 0;
 305   if (fPMD) hits = fPMD->Hits();
 306
 307   // Start loop on tracks in the hits containers
 308
 309   for (Int_t track=0; track<ntracks;track++)
 310     {
 311       gAlice->GetMCApp()->ResetHits();
 312       treeH->GetEvent(track);
 313       if (fPMD)
 314         {
 315           npmd = hits->GetEntriesFast();
 316           for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
 317             {
 318               fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
 319               trackno = fPMDHit->GetTrack();
 320               //  get kinematics of the particles
 321
 322               TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
 323               trackpid  = mparticle->GetPdgCode();
 324               Int_t  ks = mparticle->GetStatusCode();
 325               Int_t imo;
 326               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
 327
 328               if (mparticle->GetFirstMother() == -1)
 329                 {
 330                   tracknoOld  = trackno;
 331                   trackpidOld = trackpid;
 332                   statusOld   = -1;
 333                 }
 334               Int_t igstatus = 0;
 335
 336               Int_t trnotemp = trackno;    // Modified on 25th Nov 2009
 337               if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
 338                 vx = mparticle->Vx();
 339                 vy = mparticle->Vy();
 340                 vz = mparticle->Vz();
 341
 342                 if(trackpid==kGamma||trackpid==11||trackpid==-11||
 343                    trackpid==kPi0)igstatus=1;
 344               }
 345
 346
 347               while(((imo = mparticle->GetFirstMother()) >= 0)&&
 348                     (ks = mparticle->GetStatusCode() <1) )
 349                 {
 350                   mparticle =  gAlice->GetMCApp()->Particle(imo);
 351                   trackpid = mparticle->GetPdgCode();
 352                   ks = mparticle->GetStatusCode();
 353                   vx = mparticle->Vx();
 354                   vy = mparticle->Vy();
 355                   vz = mparticle->Vz();
 356
 357                   // Modified on 25th Nov 2009
 358
 359                   trnotemp = trackno;
 360                   if(trackpid == 111)
 361                     {
 362                       trackno = trnotemp;
 363                     }
 364                   if(trackpid != 111)
 365                     {
 366                       trackno=imo;
 367                     }
 368                   // end of modification on 25th Nov 2009
 369                 }
 370
 371               if(trackpid==kGamma||trackpid==11||trackpid==-11||
 372                  trackpid==kPi0)igstatus=1;
 373               mtrackpid=trackpid;
 374               mtrackno=trackno;
 375               trackpid=trackpidOld;
 376               trackno=tracknoOld;
 377
 378               //-----------------end of modification----------------
 379               xPos = fPMDHit->X();
 380               yPos = fPMDHit->Y();
 381               zPos = fPMDHit->Z();
 382
 383               edep       = fPMDHit->GetEnergy();
 384               Int_t vol1 = fPMDHit->GetVolume(1); // Column
 385               Int_t vol2 = fPMDHit->GetVolume(2); // Row
 386               Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
 387
 388
 389               // -----------------------------------------//
 390               // In new geometry after adding electronics //
 391               // For Super Module 1 & 2                   //
 392               //  nrow = 48, ncol = 96                    //
 393               // For Super Module 3 & 4                   //
 394               //  nrow = 96, ncol = 48                    //
 395               // -----------------------------------------//
 396
 397               if (vol7 < 24)
 398                 {
 399                   smnumber = vol7;
 400                 }
 401               else
 402                 {
 403                   smnumber = vol7 - 24;
 404                 }
 405               Int_t vol8 = smnumber/6 + 1;   // fake supermodule
 406
 407               if (vol8 == 1 || vol8 == 2)
 408                 {
 409                   xpad = vol2;
 410                   ypad = vol1;
 411                 }
 412               else if (vol8 == 3 || vol8 == 4)
 413                 {
 414                   xpad = vol1;
 415                   ypad = vol2;
 416                 }
 417
 418               AliDebug(2,Form("Zposition = %f Edeposition = %f",zPos,edep));
 419
 420               if (vol7 < 24)
 421                 {
 422                   // PRE
 423                   fDetNo = 0;
 424                 }
 425               else
 426                 {
 427                   // CPV
 428                   fDetNo = 1;
 429                 }
 430
 431               Int_t smn = smnumber;
 432               Int_t ixx = xpad     - 1;
 433               Int_t iyy = ypad     - 1;
 434               if (fDetNo == 0)
 435                 {
 436                   fPRE[smn][ixx][iyy] += edep;
 437                   fPRECounter[smn][ixx][iyy]++;
 438
 439                   AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 440                   fCell.Add(cell);
 441                 }
 442               else if(fDetNo == 1)
 443                 {
 444                   fCPV[smn][ixx][iyy] += edep;
 445                   fCPVCounter[smn][ixx][iyy]++;
 446                   AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 447                   fCPVCell.Add(cpvcell);
 448                 }
 449             }
 450         }
 451     } // Track Loop ended
 452   TrackAssignment2CPVCell();
 453   TrackAssignment2Cell();
 454   ResetCell();
 455
 456   Float_t deltaE      = 0.;
 457   Int_t   detno       = 0;
 458   Int_t   trno        = -1;
 459   Int_t   trpid       = -99;
 460
 461   for (Int_t idet = 0; idet < 2; idet++)
 462     {
 463       for (Int_t ism = 0; ism < fgkTotUM; ism++)
 464         {
 465           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
 466             {
 467               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
 468                 {
 469                   if (idet == 0)
 470                     {
 471                       deltaE = fPRE[ism][jrow][kcol];
 472                       trno   = fPRETrackNo[ism][jrow][kcol];
 473                       detno = 0;
 474                     }
 475                   else if (idet == 1)
 476                     {
 477                       deltaE = fCPV[ism][jrow][kcol];
 478                       trno   = fCPVTrackNo[ism][jrow][kcol];
 479                       detno = 1;
 480                     }
 481                   if (deltaE > 0.)
 482                     {
 483                       // Natasha
 484                       TParticle *mparticle = gAlice->GetMCApp()->Particle(trno);
 485                       trpid = mparticle->GetPdgCode();
 486                       AddSDigit(trno,trpid,detno,ism,jrow,kcol,deltaE);
 487                     }
 488                 }
 489             }
 490           treeS->Fill();
 491           ResetSDigit();
 492         }
 493     }
 494   fPMDLoader->WriteSDigits("OVERWRITE");
 495   ResetCellADC();
 496 }
 497 //____________________________________________________________________________
 498
 499 void AliPMDDigitizer::Hits2Digits(Int_t ievt)
 500 {
 501   // This reads the PMD Hits tree and assigns the right track number
 502   // to a cell and stores in the digits tree
 503   //
 504   const Int_t kPi0 = 111;
 505   const Int_t kGamma = 22;
 506   Int_t npmd;
 507   Int_t trackno;
 508   Int_t smnumber;
 509   Int_t trackpid;
 510   Int_t mtrackno;
 511   Int_t mtrackpid;
 512
 513   Float_t xPos, yPos, zPos;
 514   Int_t xpad = -1, ypad = -1;
 515   Float_t edep;
 516   Float_t vx = -999.0, vy = -999.0, vz = -999.0;
 517
 518   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
 519   ResetDigit();
 520
 521   AliDebug(1,Form("Event Number =  %d",ievt));
 522   Int_t nparticles = fRunLoader->GetHeader()->GetNtrack();
 523   AliDebug(1,Form("Number of Particles = %d", nparticles));
 524
 525   fRunLoader->GetEvent(ievt);
 526   // ------------------------------------------------------- //
 527   // Pointer to specific detector hits.
 528   // Get pointers to Alice detectors and Hits containers
 529
 530   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
 531   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 532
 533   if (fPMDLoader == 0x0)
 534     {
 535       AliError("Can not find PMD or PMDLoader");
 536     }
 537   TTree* treeH = fPMDLoader->TreeH();
 538   Int_t ntracks    = (Int_t) treeH->GetEntries();
 539   AliDebug(1,Form("Number of Tracks in the TreeH = %d", ntracks));
 540   fPMDLoader->LoadDigits("recreate");
 541   TTree* treeD = fPMDLoader->TreeD();
 542   if (treeD == 0x0)
 543     {
 544       fPMDLoader->MakeTree("D");
 545       treeD = fPMDLoader->TreeD();
 546     }
 547   Int_t bufsize = 16000;
 548   treeD->Branch("PMDDigit", &fDigits, bufsize);
 549
 550   TClonesArray* hits = 0;
 551   if (fPMD) hits = fPMD->Hits();
 552
 553   // Start loop on tracks in the hits containers
 554
 555   for (Int_t track=0; track<ntracks;track++)
 556     {
 557       gAlice->GetMCApp()->ResetHits();
 558       treeH->GetEvent(track);
 559
 560       if (fPMD)
 561         {
 562           npmd = hits->GetEntriesFast();
 563           for (Int_t ipmd = 0; ipmd < npmd; ipmd++)
 564             {
 565               fPMDHit = (AliPMDhit*) hits->UncheckedAt(ipmd);
 566               trackno = fPMDHit->GetTrack();
 567
 568               //  get kinematics of the particles
 569
 570               TParticle* mparticle = gAlice->GetMCApp()->Particle(trackno);
 571               trackpid  = mparticle->GetPdgCode();
 572               Int_t  ks = mparticle->GetStatusCode();
 573               Int_t imo;
 574               Int_t tracknoOld=0, trackpidOld=0, statusOld = 0;
 575               if (mparticle->GetFirstMother() == -1)
 576                 {
 577                   tracknoOld  = trackno;
 578                   trackpidOld = trackpid;
 579                   statusOld   = -1;
 580                 }
 581
 582               Int_t igstatus = 0;
 583
 584               Int_t trnotemp = trackno;   // modified on 25th Nov 2009
 585               if(ks==1||(imo = mparticle->GetFirstMother())<0 ){
 586                 vx = mparticle->Vx();
 587                 vy = mparticle->Vy();
 588                 vz = mparticle->Vz();
 589
 590                 if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
 591                   igstatus=1;
 592               }
 593
 594
 595               while(((imo = mparticle->GetFirstMother()) >= 0)&&
 596                     (ks = mparticle->GetStatusCode() <1) )
 597                 {
 598                   mparticle =  gAlice->GetMCApp()->Particle(imo);
 599                   trackpid = mparticle->GetPdgCode();
 600                   ks = mparticle->GetStatusCode();
 601                   vx = mparticle->Vx();
 602                   vy = mparticle->Vy();
 603                   vz = mparticle->Vz();
 604
 605                   // Modified on 25th Nov 2009
 606
 607                   trnotemp = trackno;
 608                   if(trackpid == 111)
 609                     {
 610                       trackno = trnotemp;
 611                     }
 612                   if(trackpid != 111)
 613                     {
 614                       trackno=imo;
 615                     }
 616                 }
 617
 618               if(trackpid==kGamma||trackpid==11||trackpid==-11||trackpid==kPi0)
 619                 igstatus=1;
 620               mtrackpid=trackpid;
 621               mtrackno=trackno;
 622               trackpid=trackpidOld;
 623               trackno=tracknoOld;
 624
 625               xPos = fPMDHit->X();
 626               yPos = fPMDHit->Y();
 627               zPos = fPMDHit->Z();
 628               edep       = fPMDHit->GetEnergy();
 629               Int_t vol1 = fPMDHit->GetVolume(1); // Column
 630               Int_t vol2 = fPMDHit->GetVolume(2); // Row
 631               Int_t vol7 = fPMDHit->GetVolume(4); // Serial Module No
 632
 633               // -----------------------------------------//
 634               // In new geometry after adding electronics //
 635               // For Super Module 1 & 2                   //
 636               //  nrow = 48, ncol = 96                    //
 637               // For Super Module 3 & 4                   //
 638               //  nrow = 96, ncol = 48                    //
 639               // -----------------------------------------//
 640
 641               if (vol7 < 24)
 642                 {
 643                   smnumber = vol7;
 644                 }
 645               else
 646                 {
 647                   smnumber = vol7 - 24;
 648                 }
 649               Int_t vol8 = smnumber/6 + 1;    // fake supermodule
 650
 651               if (vol8 == 1 || vol8 == 2)
 652                 {
 653                   xpad = vol2;
 654                   ypad = vol1;
 655                 }
 656               else if (vol8 == 3 || vol8 == 4)
 657                 {
 658                   xpad = vol1;
 659                   ypad = vol2;
 660                 }
 661
 662               AliDebug(2,Form("ZPosition = %f Edeposition = %d",zPos,edep));
 663
 664               if (vol7 < 24)
 665                 {
 666                   // PRE
 667                   fDetNo = 0;
 668                 }
 669               else
 670                 {
 671                   fDetNo = 1;
 672                 }
 673               /*
 674                 if (zPos < fZPos)
 675                 {
 676                 // CPV
 677                 fDetNo = 1;
 678                 }
 679                 else if (zPos > fZPos)
 680                 {
 681                 // PMD
 682                 fDetNo = 0;
 683                 }
 684               */
 685               Int_t smn = smnumber;
 686               Int_t ixx = xpad     - 1;
 687               Int_t iyy = ypad     - 1;
 688               if (fDetNo == 0)
 689                 {
 690                   fPRE[smn][ixx][iyy] += edep;
 691                   fPRECounter[smn][ixx][iyy]++;
 692
 693                   AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 694
 695                   fCell.Add(cell);
 696                 }
 697               else if(fDetNo == 1)
 698                 {
 699                   fCPV[smn][ixx][iyy] += edep;
 700                   fCPVCounter[smn][ixx][iyy]++;
 701                   AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 702                   fCPVCell.Add(cpvcell);
 703                 }
 704             }
 705         }
 706     } // Track Loop ended
 707   TrackAssignment2CPVCell();
 708   TrackAssignment2Cell();
 709   ResetCell();
 710
 711   Float_t gain1;
 712   Float_t adc;
 713   Float_t deltaE = 0.;
 714   Int_t detno = 0;
 715   Int_t trno = 1;
 716   Int_t trpid = -99;
 717
 718   for (Int_t idet = 0; idet < 2; idet++)
 719   {
 720       for (Int_t ism = 0; ism < fgkTotUM; ism++)
 721       {
 722           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
 723           {
 724               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
 725               {
 726                   if (idet == 0)
 727                   {
 728                       deltaE = fPRE[ism][jrow][kcol];
 729                       trno   = fPRETrackNo[ism][jrow][kcol];
 730                       detno  = 0;
 731                   }
 732                   else if (idet == 1)
 733                   {
 734                       deltaE = fCPV[ism][jrow][kcol];
 735                       trno   = fCPVTrackNo[ism][jrow][kcol];
 736                       detno  = 1;
 737                   }
 738                   if (deltaE > 0.)
 739                   {
 740                       MeV2ADC(deltaE,adc);
 741
 742                       // To decalibrate the adc values
 743                       //
 744                       gain1 = Gain(idet,ism,jrow,kcol);
 745                       if (gain1 != 0.)
 746                       {
 747                           Int_t adcDecalib = (Int_t)(adc/gain1);
 748                           adc = (Float_t) adcDecalib;
 749                       }
 750                       else if(gain1 == 0.)
 751                       {
 752                           adc = 0.;
 753                       }
 754
 755                       // Pedestal Decalibration
 756                       Int_t   pedmeanrms =
 757                           fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
 758                       Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 759                       Float_t pedrms     = (Float_t)pedrms1/10.;
 760                       Float_t pedmean    =
 761                           (Float_t) (pedmeanrms - pedrms1)/1000.0;
 762                       if (adc > 0.)
 763                       {
 764                           adc += (pedmean + 3.0*pedrms);
 765                           TParticle *mparticle
 766                             = gAlice->GetMCApp()->Particle(trno);
 767                           trpid = mparticle->GetPdgCode();
 768
 769                           AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
 770                       }
 771                   }
 772               } // column loop
 773           } // row    loop
 774           treeD->Fill();
 775           ResetDigit();
 776       } // supermodule loop
 777   } // detector loop
 778
 779   fPMDLoader->WriteDigits("OVERWRITE");
 780   ResetCellADC();
 781
 782 }
 783 //____________________________________________________________________________
 784
 785
 786 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
 787 {
 788   // This reads the PMD sdigits tree and converts energy deposition
 789   // in a cell to ADC and stores in the digits tree
 790   //
 791
 792   fRunLoader->GetEvent(ievt);
 793
 794   TTree* treeS = fPMDLoader->TreeS();
 795   AliPMDsdigit  *pmdsdigit;
 796   TBranch *branch = treeS->GetBranch("PMDSDigit");
 797   if(!branch)
 798     {
 799       AliError("PMD Sdigit branch does not exist");
 800       return;
 801     }
 802   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
 803   branch->SetAddress(&fSDigits);
 804
 805   TTree* treeD = fPMDLoader->TreeD();
 806   if (treeD == 0x0)
 807     {
 808       fPMDLoader->MakeTree("D");
 809       treeD = fPMDLoader->TreeD();
 810     }
 811   Int_t bufsize = 16000;
 812   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
 813   treeD->Branch("PMDDigit", &fDigits, bufsize);
 814
 815   Int_t   trno, trpid, det, smn;
 816   Int_t   irow, icol;
 817   Float_t edep, adc;
 818
 819   Int_t nmodules = (Int_t) treeS->GetEntries();
 820   AliDebug(1,Form("Number of modules = %d",nmodules));
 821
 822   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 823     {
 824       treeS->GetEntry(imodule);
 825       Int_t nentries = fSDigits->GetLast();
 826       AliDebug(2,Form("Number of entries per module = %d",nentries+1));
 827       for (Int_t ient = 0; ient < nentries+1; ient++)
 828         {
 829           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
 830           trno   = pmdsdigit->GetTrackNumber();
 831           trpid  = pmdsdigit->GetTrackPid();
 832           det    = pmdsdigit->GetDetector();
 833           smn    = pmdsdigit->GetSMNumber();
 834           irow   = pmdsdigit->GetRow();
 835           icol   = pmdsdigit->GetColumn();
 836           edep   = pmdsdigit->GetCellEdep();
 837
 838           MeV2ADC(edep,adc);
 839
 840           // To decalibrte the adc values
 841           //
 842           Float_t gain1 = Gain(det,smn,irow,icol);
 843           if (gain1 != 0.)
 844           {
 845             Int_t adcDecalib = (Int_t)(adc/gain1);
 846             adc = (Float_t) adcDecalib;
 847           }
 848           else if(gain1 == 0.)
 849           {
 850               adc = 0.;
 851           }
 852           // Pedestal Decalibration
 853           Int_t   pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,irow,icol);
 854           Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 855           Float_t pedrms     = (Float_t)pedrms1/10.;
 856           Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 857           if(adc > 0.)
 858           {
 859               adc += (pedmean + 3.0*pedrms);
 860               AddDigit(trno,trpid,det,smn,irow,icol,adc);
 861           }
 862
 863         }
 864       treeD->Fill();
 865       ResetDigit();
 866     }
 867   fPMDLoader->WriteDigits("OVERWRITE");
 868
 869 }
 870 //____________________________________________________________________________
 871 void AliPMDDigitizer::Exec(Option_t *option)
 872 {
 873   // Does the event merging and digitization
 874   const char *cdeb = strstr(option,"deb");
 875   if(cdeb)
 876     {
 877       AliDebug(100," *** PMD Exec is called ***");
 878     }
 879
 880   Int_t ninputs = fManager->GetNinputs();
 881   AliDebug(1,Form("Number of files to be processed = %d",ninputs));
 882   ResetCellADC();
 883
 884   for (Int_t i = 0; i < ninputs; i++)
 885     {
 886       Int_t troffset = fManager->GetMask(i);
 887       MergeSDigits(i, troffset);
 888     }
 889
 890   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 891   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
 892   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 893   if (fPMDLoader == 0x0)
 894     {
 895       AliError("Can not find PMD or PMDLoader");
 896     }
 897   fPMDLoader->LoadDigits("update");
 898   TTree* treeD = fPMDLoader->TreeD();
 899   if (treeD == 0x0)
 900     {
 901       fPMDLoader->MakeTree("D");
 902       treeD = fPMDLoader->TreeD();
 903     }
 904   Int_t bufsize = 16000;
 905   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
 906   treeD->Branch("PMDDigit", &fDigits, bufsize);
 907
 908   Float_t adc;
 909   Float_t deltaE = 0.;
 910   Int_t detno    = 0;
 911   Int_t trno     = 1;
 912   Int_t trpid    = -99;
 913
 914   for (Int_t idet = 0; idet < 2; idet++)
 915     {
 916       for (Int_t ism = 0; ism < fgkTotUM; ism++)
 917         {
 918           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
 919             {
 920               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
 921                 {
 922                   if (idet == 0)
 923                     {
 924                       deltaE = fPRE[ism][jrow][kcol];
 925                       trno   = fPRETrackNo[ism][jrow][kcol];
 926                       trpid  = fPRETrackPid[ism][jrow][kcol];
 927                       detno  = 0;
 928                     }
 929                   else if (idet == 1)
 930                     {
 931                       deltaE = fCPV[ism][jrow][kcol];
 932                       trno   = fCPVTrackNo[ism][jrow][kcol];
 933                       trpid  = fCPVTrackPid[ism][jrow][kcol];
 934                       detno  = 1;
 935                     }
 936                   if (deltaE > 0.)
 937                     {
 938                       MeV2ADC(deltaE,adc);
 939
 940                       //
 941                       // Gain decalibration
 942                       //
 943                       Float_t gain1 = Gain(idet,ism,jrow,kcol);
 944
 945                       if (gain1 != 0.)
 946                       {
 947                           Int_t adcDecalib = (Int_t)(adc/gain1);
 948                           adc = (Float_t) adcDecalib;
 949                       }
 950                       else if(gain1 == 0.)
 951                       {
 952                           adc = 0.;
 953                       }
 954                       // Pedestal Decalibration
 955                       Int_t   pedmeanrms =
 956                           fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
 957                       Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 958                       Float_t pedrms     = (Float_t)pedrms1/10.;
 959                       Float_t pedmean    =
 960                           (Float_t) (pedmeanrms - pedrms1)/1000.0;
 961                       if (adc > 0.)
 962                       {
 963                           adc += (pedmean + 3.0*pedrms);
 964                           AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
 965                       }
 966
 967                     }
 968                 } // column loop
 969             } // row    loop
 970           treeD->Fill();
 971           ResetDigit();
 972         } // supermodule loop
 973     } // detector loop
 974   fPMDLoader->WriteDigits("OVERWRITE");
 975   fPMDLoader->UnloadDigits();
 976   ResetCellADC();
 977 }
 978 //____________________________________________________________________________
 979 void AliPMDDigitizer::TrackAssignment2CPVCell()
 980 {
 981   // This block assigns the cell id when there are
 982   // multiple tracks in a cell according to the
 983   // energy deposition
 984   // This method added by Ajay
 985   Bool_t jsort = false;
 986
 987   Int_t i, j, k;
 988
 989   Int_t   *status1;
 990   Int_t   *status2;
 991   Int_t   *trnarray;
 992   Float_t *fracEdp;
 993   Float_t *trEdp;
 994
 995   Int_t   ****cpvTrack;
 996   Float_t ****cpvEdep;
 997
 998   cpvTrack = new Int_t ***[fgkTotUM];
 999   cpvEdep  = new Float_t ***[fgkTotUM];
1000   for (i=0; i<fgkTotUM; i++)
1001     {
1002       cpvTrack[i] = new Int_t **[fgkRow];
1003       cpvEdep[i]  = new Float_t **[fgkRow];
1004     }
1005
1006   for (i = 0; i < fgkTotUM; i++)
1007     {
1008       for (j = 0; j < fgkRow; j++)
1009         {
1010           cpvTrack[i][j] = new Int_t *[fgkCol];
1011           cpvEdep[i][j]  = new Float_t *[fgkCol];
1012         }
1013     }
1014   for (i = 0; i < fgkTotUM; i++)
1015     {
1016       for (j = 0; j < fgkRow; j++)
1017         {
1018           for (k = 0; k < fgkCol; k++)
1019             {
1020               Int_t nn = fCPVCounter[i][j][k];
1021               if(nn > 0)
1022                 {
1023                   cpvTrack[i][j][k] = new Int_t[nn];
1024                   cpvEdep[i][j][k] = new Float_t[nn];
1025                 }
1026               else
1027                 {
1028                   nn = 1;
1029                   cpvTrack[i][j][k] = new Int_t[nn];
1030                   cpvEdep[i][j][k] = new Float_t[nn];
1031                 }
1032               fCPVCounter[i][j][k] = 0;
1033             }
1034         }
1035     }
1036
1037
1038   Int_t nentries = fCPVCell.GetEntries();
1039
1040   Int_t   mtrackno, ism, ixp, iyp;
1041   Float_t edep;
1042   for (i = 0; i < nentries; i++)
1043     {
1044       AliPMDcell* cpvcell = (AliPMDcell*)fCPVCell.UncheckedAt(i);
1045
1046       mtrackno = cpvcell->GetTrackNumber();
1047       ism      = cpvcell->GetSMNumber();
1048       ixp      = cpvcell->GetX();
1049       iyp      = cpvcell->GetY();
1050       edep     = cpvcell->GetEdep();
1051       Int_t nn = fCPVCounter[ism][ixp][iyp];
1052       cpvTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1053       cpvEdep[ism][ixp][iyp][nn] = edep;
1054       fCPVCounter[ism][ixp][iyp]++;
1055     }
1056
1057   Int_t iz, il;
1058   Int_t im, ix, iy;
1059   Int_t nn;
1060   for (im=0; im<fgkTotUM; im++)
1061     {
1062       for (ix=0; ix<fgkRow; ix++)
1063         {
1064           for (iy=0; iy<fgkCol; iy++)
1065             {
1066               nn = fCPVCounter[im][ix][iy];
1067               if (nn > 1)
1068                 {
1069                   // This block handles if a cell is fired
1070                   // many times by many tracks
1071                   status1  = new Int_t[nn];
1072                   status2  = new Int_t[nn];
1073                   trnarray = new Int_t[nn];
1074                   for (iz = 0; iz < nn; iz++)
1075                     {
1076                       status1[iz] = cpvTrack[im][ix][iy][iz];
1077                     }
1078                   TMath::Sort(nn,status1,status2,jsort);
1079                   Int_t trackOld = -99999;
1080                   Int_t track, trCount = 0;
1081                   for (iz = 0; iz < nn; iz++)
1082                     {
1083                       track = status1[status2[iz]];
1084                       if (trackOld != track)
1085                         {
1086                           trnarray[trCount] = track;
1087                           trCount++;
1088                         }
1089                       trackOld = track;
1090                     }
1091                   delete [] status1;
1092                   delete [] status2;
1093                   Float_t totEdp = 0.;
1094                   trEdp = new Float_t[trCount];
1095                   fracEdp = new Float_t[trCount];
1096                   for (il = 0; il < trCount; il++)
1097                     {
1098                       trEdp[il] = 0.;
1099                       track = trnarray[il];
1100                       for (iz = 0; iz < nn; iz++)
1101                         {
1102                           if (track == cpvTrack[im][ix][iy][iz])
1103                             {
1104                               trEdp[il] += cpvEdep[im][ix][iy][iz];
1105                             }
1106                         }
1107                       totEdp += trEdp[il];
1108                     }
1109                   Int_t ilOld = 0;
1110                   Float_t fracOld = 0.;
1111
1112                   for (il = 0; il < trCount; il++)
1113                     {
1114                       fracEdp[il] = trEdp[il]/totEdp;
1115                       if (fracOld < fracEdp[il])
1116                         {
1117                           fracOld = fracEdp[il];
1118                           ilOld = il;
1119                         }
1120                     }
1121                   fCPVTrackNo[im][ix][iy] = trnarray[ilOld];
1122                   delete [] fracEdp;
1123                   delete [] trEdp;
1124                   delete [] trnarray;
1125                 }
1126               else if (nn == 1)
1127                 {
1128                   // This only handles if a cell is fired
1129                   // by only one track
1130
1131                   fCPVTrackNo[im][ix][iy] = cpvTrack[im][ix][iy][0];
1132
1133                 }
1134               else if (nn ==0)
1135                 {
1136                   // This is if no cell is fired
1137                   fCPVTrackNo[im][ix][iy] = -999;
1138                 }
1139             } // end of iy
1140         } // end of ix
1141     } // end of im
1142
1143   // Delete all the pointers
1144
1145  for (i = 0; i < fgkTotUM; i++)
1146     {
1147       for (j = 0; j < fgkRow; j++)
1148         {
1149           for (k = 0; k < fgkCol; k++)
1150             {
1151               delete []cpvTrack[i][j][k];
1152               delete []cpvEdep[i][j][k];
1153             }
1154         }
1155     }
1156
1157   for (i = 0; i < fgkTotUM; i++)
1158     {
1159       for (j = 0; j < fgkRow; j++)
1160         {
1161           delete [] cpvTrack[i][j];
1162           delete [] cpvEdep[i][j];
1163         }
1164     }
1165
1166   for (i = 0; i < fgkTotUM; i++)
1167     {
1168       delete [] cpvTrack[i];
1169       delete [] cpvEdep[i];
1170     }
1171   delete [] cpvTrack;
1172   delete [] cpvEdep;
1173
1174   //
1175   // End of the cell id assignment
1176   //
1177 }
1178 //____________________________________________________________________________
1179
1180 void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
1181 {
1182   // merging sdigits
1183   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(filenumber));
1184   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
1185   fPMDLoader->LoadSDigits("read");
1186   TTree* treeS = fPMDLoader->TreeS();
1187   AliPMDsdigit  *pmdsdigit;
1188   TBranch *branch = treeS->GetBranch("PMDSDigit");
1189   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1190   branch->SetAddress(&fSDigits);
1191
1192   Int_t   itrackno, itrackpid, idet, ism;
1193   Int_t   ixp, iyp;
1194   Float_t edep;
1195   Int_t nmodules = (Int_t) treeS->GetEntries();
1196   AliDebug(1,Form("Number of Modules in the treeS = %d",nmodules));
1197   AliDebug(1,Form("Track Offset = %d",troffset));
1198   for (Int_t imodule = 0; imodule < nmodules; imodule++)
1199     {
1200       treeS->GetEntry(imodule);
1201       Int_t nentries = fSDigits->GetLast();
1202       AliDebug(2,Form("Number of Entries per Module = %d",nentries));
1203       for (Int_t ient = 0; ient < nentries+1; ient++)
1204         {
1205           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
1206           itrackno  = pmdsdigit->GetTrackNumber();
1207           itrackpid = pmdsdigit->GetTrackPid();
1208           idet      = pmdsdigit->GetDetector();
1209           ism       = pmdsdigit->GetSMNumber();
1210           ixp       = pmdsdigit->GetRow();
1211           iyp       = pmdsdigit->GetColumn();
1212           edep      = pmdsdigit->GetCellEdep();
1213           if (idet == 0)
1214             {
1215               if (fPRE[ism][ixp][iyp] < edep)
1216                 {
1217                   fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
1218                   fPRETrackPid[ism][ixp][iyp] = itrackpid;
1219                 }
1220               fPRE[ism][ixp][iyp] += edep;
1221             }
1222           else if (idet == 1)
1223             {
1224               if (fCPV[ism][ixp][iyp] < edep)
1225                 {
1226                   fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
1227                   fCPVTrackPid[ism][ixp][iyp] = itrackpid;
1228                 }
1229               fCPV[ism][ixp][iyp] += edep;
1230             }
1231         }
1232     }
1233
1234 }
1235 // ----------------------------------------------------------------------
1236 void AliPMDDigitizer::TrackAssignment2Cell()
1237 {
1238   //
1239   // This block assigns the cell id when there are
1240   // multiple tracks in a cell according to the
1241   // energy deposition
1242   //
1243   Bool_t jsort = false;
1244
1245   Int_t i, j, k;
1246
1247   Int_t   *status1;
1248   Int_t   *status2;
1249   Int_t   *trnarray;
1250   Float_t *fracEdp;
1251   Float_t *trEdp;
1252
1253   Int_t   ****pmdTrack;
1254   Float_t ****pmdEdep;
1255
1256   pmdTrack = new Int_t ***[fgkTotUM];
1257   pmdEdep  = new Float_t ***[fgkTotUM];
1258   for (i=0; i<fgkTotUM; i++)
1259     {
1260       pmdTrack[i] = new Int_t **[fgkRow];
1261       pmdEdep[i]  = new Float_t **[fgkRow];
1262     }
1263
1264   for (i = 0; i < fgkTotUM; i++)
1265     {
1266       for (j = 0; j < fgkRow; j++)
1267         {
1268           pmdTrack[i][j] = new Int_t *[fgkCol];
1269           pmdEdep[i][j]  = new Float_t *[fgkCol];
1270         }
1271     }
1272
1273   for (i = 0; i < fgkTotUM; i++)
1274     {
1275       for (j = 0; j < fgkRow; j++)
1276         {
1277           for (k = 0; k < fgkCol; k++)
1278             {
1279               Int_t nn = fPRECounter[i][j][k];
1280               if(nn > 0)
1281                 {
1282                   pmdTrack[i][j][k] = new Int_t[nn];
1283                   pmdEdep[i][j][k] = new Float_t[nn];
1284                 }
1285               else
1286                 {
1287                   nn = 1;
1288                   pmdTrack[i][j][k] = new Int_t[nn];
1289                   pmdEdep[i][j][k] = new Float_t[nn];
1290                 }
1291               fPRECounter[i][j][k] = 0;
1292             }
1293         }
1294     }
1295
1296
1297   Int_t nentries = fCell.GetEntries();
1298
1299   Int_t   mtrackno, ism, ixp, iyp;
1300   Float_t edep;
1301
1302   for (i = 0; i < nentries; i++)
1303     {
1304       AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
1305
1306       mtrackno = cell->GetTrackNumber();
1307       ism      = cell->GetSMNumber();
1308       ixp      = cell->GetX();
1309       iyp      = cell->GetY();
1310       edep     = cell->GetEdep();
1311       Int_t nn = fPRECounter[ism][ixp][iyp];
1312       pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1313       pmdEdep[ism][ixp][iyp][nn] = edep;
1314       fPRECounter[ism][ixp][iyp]++;
1315     }
1316
1317   Int_t iz, il;
1318   Int_t im, ix, iy;
1319   Int_t nn;
1320
1321   for (im=0; im<fgkTotUM; im++)
1322     {
1323       for (ix=0; ix<fgkRow; ix++)
1324         {
1325           for (iy=0; iy<fgkCol; iy++)
1326             {
1327               nn = fPRECounter[im][ix][iy];
1328               if (nn > 1)
1329                 {
1330                   // This block handles if a cell is fired
1331                   // many times by many tracks
1332                   status1  = new Int_t[nn];
1333                   status2  = new Int_t[nn];
1334                   trnarray = new Int_t[nn];
1335                   for (iz = 0; iz < nn; iz++)
1336                     {
1337                       status1[iz] = pmdTrack[im][ix][iy][iz];
1338                     }
1339                   TMath::Sort(nn,status1,status2,jsort);
1340                   Int_t trackOld = -99999;
1341                   Int_t track, trCount = 0;
1342                   for (iz = 0; iz < nn; iz++)
1343                     {
1344                       track = status1[status2[iz]];
1345                       if (trackOld != track)
1346                         {
1347                           trnarray[trCount] = track;
1348                           trCount++;
1349                         }
1350                       trackOld = track;
1351                     }
1352                   delete [] status1;
1353                   delete [] status2;
1354                   Float_t totEdp = 0.;
1355                   trEdp = new Float_t[trCount];
1356                   fracEdp = new Float_t[trCount];
1357                   for (il = 0; il < trCount; il++)
1358                     {
1359                       trEdp[il] = 0.;
1360                       track = trnarray[il];
1361                       for (iz = 0; iz < nn; iz++)
1362                         {
1363                           if (track == pmdTrack[im][ix][iy][iz])
1364                             {
1365                               trEdp[il] += pmdEdep[im][ix][iy][iz];
1366                             }
1367                         }
1368                       totEdp += trEdp[il];
1369                     }
1370                   Int_t ilOld = 0;
1371                   Float_t fracOld = 0.;
1372
1373                   for (il = 0; il < trCount; il++)
1374                     {
1375                       fracEdp[il] = trEdp[il]/totEdp;
1376                       if (fracOld < fracEdp[il])
1377                         {
1378                           fracOld = fracEdp[il];
1379                           ilOld = il;
1380                         }
1381                     }
1382                   fPRETrackNo[im][ix][iy] = trnarray[ilOld];
1383                   delete [] fracEdp;
1384                   delete [] trEdp;
1385                   delete [] trnarray;
1386                 }
1387               else if (nn == 1)
1388                 {
1389                   // This only handles if a cell is fired
1390                   // by only one track
1391
1392                   fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
1393
1394                 }
1395               else if (nn ==0)
1396                 {
1397                   // This is if no cell is fired
1398                   fPRETrackNo[im][ix][iy] = -999;
1399                 }
1400             } // end of iy
1401         } // end of ix
1402     } // end of im
1403
1404   // Delete all the pointers
1405
1406   for (i = 0; i < fgkTotUM; i++)
1407     {
1408       for (j = 0; j < fgkRow; j++)
1409         {
1410           for (k = 0; k < fgkCol; k++)
1411             {
1412               delete [] pmdTrack[i][j][k];
1413               delete [] pmdEdep[i][j][k];
1414             }
1415         }
1416     }
1417
1418   for (i = 0; i < fgkTotUM; i++)
1419     {
1420       for (j = 0; j < fgkRow; j++)
1421         {
1422           delete [] pmdTrack[i][j];
1423           delete [] pmdEdep[i][j];
1424         }
1425     }
1426
1427   for (i = 0; i < fgkTotUM; i++)
1428     {
1429       delete [] pmdTrack[i];
1430       delete [] pmdEdep[i];
1431     }
1432   delete [] pmdTrack;
1433   delete [] pmdEdep;
1434   //
1435   // End of the cell id assignment
1436   //
1437 }
1438 //____________________________________________________________________________
1439 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
1440 {
1441   // This converts the simulated edep to ADC according to the
1442   // Test Beam Data
1443   //PS Test in May 2009, Voltage @ 1350 V
1444   // KeV - ADC conversion for 12bit ADC
1445   // MPV data used for the fit and taken here
1446
1447   const Float_t kConstant   = -0.1602;
1448   const Float_t kErConstant =  0.9914;
1449   const Float_t kSlope      = 77.47;
1450   const Float_t kErSlope    =  3.16;
1451
1452   Float_t cons   = gRandom->Gaus(kConstant,kErConstant);
1453   Float_t slop   = gRandom->Gaus(kSlope,kErSlope);
1454
1455   Float_t adc12bit = slop*mev*0.001 + cons;
1456
1457   if (adc12bit < 0.) adc12bit = 0.;
1458
1459   if(adc12bit < 1600.0)
1460     {
1461       adc = (Float_t) adc12bit;
1462     }
1463   else if (adc12bit >= 1600.0)
1464     {
1465       adc = 1600.0;
1466     }
1467 }
1468 //____________________________________________________________________________
1469 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t trpid, Int_t det,
1470                                 Int_t smnumber, Int_t irow, Int_t icol,
1471                                 Float_t adc)
1472 {
1473   // Add SDigit
1474   //
1475   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1476   TClonesArray &lsdigits = *fSDigits;
1477   new(lsdigits[fNsdigit++])  AliPMDsdigit(trnumber,trpid,det,smnumber,irow,icol,adc);
1478 }
1479 //____________________________________________________________________________
1480
1481 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t trpid, Int_t det,
1482                                Int_t smnumber, Int_t irow, Int_t icol,
1483                                Float_t adc)
1484 {
1485   // Add Digit
1486   //
1487   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
1488   TClonesArray &ldigits = *fDigits;
1489   new(ldigits[fNdigit++]) AliPMDdigit(trnumber,trpid, det,smnumber,irow,icol,adc);
1490 }
1491 //____________________________________________________________________________
1492
1493 void AliPMDDigitizer::SetZPosition(Float_t zpos)
1494 {
1495   fZPos = zpos;
1496 }
1497 //____________________________________________________________________________
1498 Float_t AliPMDDigitizer::GetZPosition() const
1499 {
1500   return fZPos;
1501 }
1502 //____________________________________________________________________________
1503
1504 void AliPMDDigitizer::ResetCell()
1505 {
1506   // clears the cell array and also the counter
1507   //  for each cell
1508   //
1509   fCPVCell.Delete();
1510   fCell.Delete();
1511   for (Int_t i = 0; i < fgkTotUM; i++)
1512     {
1513       for (Int_t j = 0; j < fgkRow; j++)
1514         {
1515           for (Int_t k = 0; k < fgkCol; k++)
1516             {
1517               fCPVCounter[i][j][k] = 0;
1518               fPRECounter[i][j][k] = 0;
1519             }
1520         }
1521     }
1522 }
1523 //____________________________________________________________________________
1524 void AliPMDDigitizer::ResetSDigit()
1525 {
1526   // Clears SDigits
1527   fNsdigit = 0;
1528   if (fSDigits) fSDigits->Delete();
1529 }
1530 //____________________________________________________________________________
1531 void AliPMDDigitizer::ResetDigit()
1532 {
1533   // Clears Digits
1534   fNdigit = 0;
1535   if (fDigits) fDigits->Delete();
1536 }
1537 //____________________________________________________________________________
1538
1539 void AliPMDDigitizer::ResetCellADC()
1540 {
1541   // Clears individual cells edep and track number
1542   for (Int_t i = 0; i < fgkTotUM; i++)
1543     {
1544       for (Int_t j = 0; j < fgkRow; j++)
1545         {
1546           for (Int_t k = 0; k < fgkCol; k++)
1547             {
1548               fCPV[i][j][k]         = 0.;
1549               fPRE[i][j][k]         = 0.;
1550               fCPVTrackNo[i][j][k]  = 0;
1551               fPRETrackNo[i][j][k]  = 0;
1552               fCPVTrackPid[i][j][k] = -1;
1553               fPRETrackPid[i][j][k] = -1;
1554             }
1555         }
1556     }
1557 }
1558 //____________________________________________________________________________
1559
1560 void AliPMDDigitizer::UnLoad(Option_t *option)
1561 {
1562   // Unloads all the root files
1563   //
1564   const char *cS = strstr(option,"S");
1565   const char *cD = strstr(option,"D");
1566
1567   fRunLoader->UnloadgAlice();
1568   fRunLoader->UnloadHeader();
1569   fRunLoader->UnloadKinematics();
1570
1571   if (cS)
1572     {
1573       fPMDLoader->UnloadHits();
1574     }
1575   if (cD)
1576     {
1577       fPMDLoader->UnloadHits();
1578       fPMDLoader->UnloadSDigits();
1579     }
1580 }
1581
1582 //----------------------------------------------------------------------
1583 Float_t AliPMDDigitizer::Gain(Int_t det, Int_t smn, Int_t row, Int_t col) const
1584 {
1585   // returns of the gain of the cell
1586   // Added this method by ZA
1587
1588   //cout<<" I am here in gain "<<fCalibData<< "smn,row, col "<<smn
1589   //<<" "<<row<<" "<<col<<endl;
1590
1591   if(!fCalibGain) {
1592     AliError("No calibration data loaded from CDB!!!");
1593     return 1;
1594   }
1595
1596   Float_t GainFact;
1597   GainFact = fCalibGain->GetGainFact(det,smn,row,col);
1598   return GainFact;
1599 }
1600 //----------------------------------------------------------------------
1601 AliPMDCalibData* AliPMDDigitizer::GetCalibGain() const
1602 {
1603   // The run number will be centralized in AliCDBManager,
1604   // you don't need to set it here!
1605   // Added this method by ZA
1606   // Cleaned up by Alberto
1607   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Gain");
1608
1609   if(!entry) AliFatal("Calibration object retrieval failed!");
1610
1611   AliPMDCalibData *calibdata=0;
1612   if (entry) calibdata = (AliPMDCalibData*) entry->GetObject();
1613
1614   if (!calibdata)  AliFatal("No calibration data from calibration database !");
1615
1616   return calibdata;
1617 }
1618 //----------------------------------------------------------------------
1619 AliPMDPedestal* AliPMDDigitizer::GetCalibPed() const
1620 {
1621   // The run number will be centralized in AliCDBManager,
1622   // you don't need to set it here!
1623
1624   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Ped");
1625
1626   if(!entry) AliFatal("Pedestal object retrieval failed!");
1627
1628   AliPMDPedestal *pedestal=0;
1629   if (entry) pedestal = (AliPMDPedestal*) entry->GetObject();
1630
1631   if (!pedestal)  AliFatal("No pedestal data from calibration database !");
1632
1633   return pedestal;
1634 }