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               //Float_t zposition = TMath::Abs(zPos);
 420               if (zPos < fZPos)
 421                 {
 422                   // CPV
 423                   fDetNo = 1;
 424                 }
 425               else if (zPos > fZPos)
 426                 {
 427                   // PMD
 428                   fDetNo = 0;
 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 (zPos < fZPos)
 665                 {
 666                   // CPV
 667                   fDetNo = 1;
 668                 }
 669               else if (zPos > fZPos)
 670                 {
 671                   // PMD
 672                   fDetNo = 0;
 673                 }
 674
 675               Int_t smn = smnumber;
 676               Int_t ixx = xpad     - 1;
 677               Int_t iyy = ypad     - 1;
 678               if (fDetNo == 0)
 679                 {
 680                   fPRE[smn][ixx][iyy] += edep;
 681                   fPRECounter[smn][ixx][iyy]++;
 682
 683                   AliPMDcell* cell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 684
 685                   fCell.Add(cell);
 686                 }
 687               else if(fDetNo == 1)
 688                 {
 689                   fCPV[smn][ixx][iyy] += edep;
 690                   fCPVCounter[smn][ixx][iyy]++;
 691                   AliPMDcell* cpvcell = new AliPMDcell(mtrackno,smn,ixx,iyy,edep);
 692                   fCPVCell.Add(cpvcell);
 693                 }
 694             }
 695         }
 696     } // Track Loop ended
 697   TrackAssignment2CPVCell();
 698   TrackAssignment2Cell();
 699   ResetCell();
 700
 701   Float_t gain1;
 702   Float_t adc;
 703   Float_t deltaE = 0.;
 704   Int_t detno = 0;
 705   Int_t trno = 1;
 706   Int_t trpid = -99;
 707
 708   for (Int_t idet = 0; idet < 2; idet++)
 709   {
 710       for (Int_t ism = 0; ism < fgkTotUM; ism++)
 711       {
 712           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
 713           {
 714               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
 715               {
 716                   if (idet == 0)
 717                   {
 718                       deltaE = fPRE[ism][jrow][kcol];
 719                       trno   = fPRETrackNo[ism][jrow][kcol];
 720                       detno  = 0;
 721                   }
 722                   else if (idet == 1)
 723                   {
 724                       deltaE = fCPV[ism][jrow][kcol];
 725                       trno   = fCPVTrackNo[ism][jrow][kcol];
 726                       detno  = 1;
 727                   }
 728                   if (deltaE > 0.)
 729                   {
 730                       MeV2ADC(deltaE,adc);
 731
 732                       // To decalibrate the adc values
 733                       //
 734                       gain1 = Gain(idet,ism,jrow,kcol);
 735                       if (gain1 != 0.)
 736                       {
 737                           Int_t adcDecalib = (Int_t)(adc/gain1);
 738                           adc = (Float_t) adcDecalib;
 739                       }
 740                       else if(gain1 == 0.)
 741                       {
 742                           adc = 0.;
 743                       }
 744
 745                       // Pedestal Decalibration
 746                       Int_t   pedmeanrms =
 747                           fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
 748                       Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 749                       Float_t pedrms     = (Float_t)pedrms1/10.;
 750                       Float_t pedmean    =
 751                           (Float_t) (pedmeanrms - pedrms1)/1000.0;
 752                       if (adc > 0.)
 753                       {
 754                           adc += (pedmean + 3.0*pedrms);
 755                           TParticle *mparticle
 756                             = gAlice->GetMCApp()->Particle(trno);
 757                           trpid = mparticle->GetPdgCode();
 758
 759                           AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
 760                       }
 761                   }
 762               } // column loop
 763           } // row    loop
 764           treeD->Fill();
 765           ResetDigit();
 766       } // supermodule loop
 767   } // detector loop
 768
 769   fPMDLoader->WriteDigits("OVERWRITE");
 770   ResetCellADC();
 771
 772 }
 773 //____________________________________________________________________________
 774
 775
 776 void AliPMDDigitizer::SDigits2Digits(Int_t ievt)
 777 {
 778   // This reads the PMD sdigits tree and converts energy deposition
 779   // in a cell to ADC and stores in the digits tree
 780   //
 781
 782   fRunLoader->GetEvent(ievt);
 783
 784   TTree* treeS = fPMDLoader->TreeS();
 785   AliPMDsdigit  *pmdsdigit;
 786   TBranch *branch = treeS->GetBranch("PMDSDigit");
 787   if(!branch)
 788     {
 789       AliError("PMD Sdigit branch does not exist");
 790       return;
 791     }
 792   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
 793   branch->SetAddress(&fSDigits);
 794
 795   TTree* treeD = fPMDLoader->TreeD();
 796   if (treeD == 0x0)
 797     {
 798       fPMDLoader->MakeTree("D");
 799       treeD = fPMDLoader->TreeD();
 800     }
 801   Int_t bufsize = 16000;
 802   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
 803   treeD->Branch("PMDDigit", &fDigits, bufsize);
 804
 805   Int_t   trno, trpid, det, smn;
 806   Int_t   irow, icol;
 807   Float_t edep, adc;
 808
 809   Int_t nmodules = (Int_t) treeS->GetEntries();
 810   AliDebug(1,Form("Number of modules = %d",nmodules));
 811
 812   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 813     {
 814       treeS->GetEntry(imodule);
 815       Int_t nentries = fSDigits->GetLast();
 816       AliDebug(2,Form("Number of entries per module = %d",nentries+1));
 817       for (Int_t ient = 0; ient < nentries+1; ient++)
 818         {
 819           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
 820           trno   = pmdsdigit->GetTrackNumber();
 821           trpid  = pmdsdigit->GetTrackPid();
 822           det    = pmdsdigit->GetDetector();
 823           smn    = pmdsdigit->GetSMNumber();
 824           irow   = pmdsdigit->GetRow();
 825           icol   = pmdsdigit->GetColumn();
 826           edep   = pmdsdigit->GetCellEdep();
 827
 828           MeV2ADC(edep,adc);
 829
 830           // To decalibrte the adc values
 831           //
 832           Float_t gain1 = Gain(det,smn,irow,icol);
 833           if (gain1 != 0.)
 834           {
 835             Int_t adcDecalib = (Int_t)(adc/gain1);
 836             adc = (Float_t) adcDecalib;
 837           }
 838           else if(gain1 == 0.)
 839           {
 840               adc = 0.;
 841           }
 842           // Pedestal Decalibration
 843           Int_t   pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,irow,icol);
 844           Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 845           Float_t pedrms     = (Float_t)pedrms1/10.;
 846           Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 847           if(adc > 0.)
 848           {
 849               adc += (pedmean + 3.0*pedrms);
 850               AddDigit(trno,trpid,det,smn,irow,icol,adc);
 851           }
 852
 853         }
 854       treeD->Fill();
 855       ResetDigit();
 856     }
 857   fPMDLoader->WriteDigits("OVERWRITE");
 858
 859 }
 860 //____________________________________________________________________________
 861 void AliPMDDigitizer::Exec(Option_t *option)
 862 {
 863   // Does the event merging and digitization
 864   const char *cdeb = strstr(option,"deb");
 865   if(cdeb)
 866     {
 867       AliDebug(100," *** PMD Exec is called ***");
 868     }
 869
 870   Int_t ninputs = fManager->GetNinputs();
 871   AliDebug(1,Form("Number of files to be processed = %d",ninputs));
 872   ResetCellADC();
 873
 874   for (Int_t i = 0; i < ninputs; i++)
 875     {
 876       Int_t troffset = fManager->GetMask(i);
 877       MergeSDigits(i, troffset);
 878     }
 879
 880   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetOutputFolderName());
 881   fPMD  = (AliPMD*)gAlice->GetDetector("PMD");
 882   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 883   if (fPMDLoader == 0x0)
 884     {
 885       AliError("Can not find PMD or PMDLoader");
 886     }
 887   fPMDLoader->LoadDigits("update");
 888   TTree* treeD = fPMDLoader->TreeD();
 889   if (treeD == 0x0)
 890     {
 891       fPMDLoader->MakeTree("D");
 892       treeD = fPMDLoader->TreeD();
 893     }
 894   Int_t bufsize = 16000;
 895   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
 896   treeD->Branch("PMDDigit", &fDigits, bufsize);
 897
 898   Float_t adc;
 899   Float_t deltaE = 0.;
 900   Int_t detno    = 0;
 901   Int_t trno     = 1;
 902   Int_t trpid    = -99;
 903
 904   for (Int_t idet = 0; idet < 2; idet++)
 905     {
 906       for (Int_t ism = 0; ism < fgkTotUM; ism++)
 907         {
 908           for (Int_t jrow = 0; jrow < fgkRow; jrow++)
 909             {
 910               for (Int_t kcol = 0; kcol < fgkCol; kcol++)
 911                 {
 912                   if (idet == 0)
 913                     {
 914                       deltaE = fPRE[ism][jrow][kcol];
 915                       trno   = fPRETrackNo[ism][jrow][kcol];
 916                       trpid  = fPRETrackPid[ism][jrow][kcol];
 917                       detno  = 0;
 918                     }
 919                   else if (idet == 1)
 920                     {
 921                       deltaE = fCPV[ism][jrow][kcol];
 922                       trno   = fCPVTrackNo[ism][jrow][kcol];
 923                       trpid  = fCPVTrackPid[ism][jrow][kcol];
 924                       detno  = 1;
 925                     }
 926                   if (deltaE > 0.)
 927                     {
 928                       MeV2ADC(deltaE,adc);
 929
 930                       //
 931                       // Gain decalibration
 932                       //
 933                       Float_t gain1 = Gain(idet,ism,jrow,kcol);
 934
 935                       if (gain1 != 0.)
 936                       {
 937                           Int_t adcDecalib = (Int_t)(adc/gain1);
 938                           adc = (Float_t) adcDecalib;
 939                       }
 940                       else if(gain1 == 0.)
 941                       {
 942                           adc = 0.;
 943                       }
 944                       // Pedestal Decalibration
 945                       Int_t   pedmeanrms =
 946                           fCalibPed->GetPedMeanRms(idet,ism,jrow,kcol);
 947                       Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 948                       Float_t pedrms     = (Float_t)pedrms1/10.;
 949                       Float_t pedmean    =
 950                           (Float_t) (pedmeanrms - pedrms1)/1000.0;
 951                       if (adc > 0.)
 952                       {
 953                           adc += (pedmean + 3.0*pedrms);
 954                           AddDigit(trno,trpid,detno,ism,jrow,kcol,adc);
 955                       }
 956
 957                     }
 958                 } // column loop
 959             } // row    loop
 960           treeD->Fill();
 961           ResetDigit();
 962         } // supermodule loop
 963     } // detector loop
 964   fPMDLoader->WriteDigits("OVERWRITE");
 965   fPMDLoader->UnloadDigits();
 966   ResetCellADC();
 967 }
 968 //____________________________________________________________________________
 969 void AliPMDDigitizer::TrackAssignment2CPVCell()
 970 {
 971   // This block assigns the cell id when there are
 972   // multiple tracks in a cell according to the
 973   // energy deposition
 974   // This method added by Ajay
 975   Bool_t jsort = false;
 976
 977   Int_t i, j, k;
 978
 979   Int_t   *status1;
 980   Int_t   *status2;
 981   Int_t   *trnarray;
 982   Float_t *fracEdp;
 983   Float_t *trEdp;
 984
 985   Int_t   ****cpvTrack;
 986   Float_t ****cpvEdep;
 987
 988   cpvTrack = new Int_t ***[fgkTotUM];
 989   cpvEdep  = new Float_t ***[fgkTotUM];
 990   for (i=0; i<fgkTotUM; i++)
 991     {
 992       cpvTrack[i] = new Int_t **[fgkRow];
 993       cpvEdep[i]  = new Float_t **[fgkRow];
 994     }
 995
 996   for (i = 0; i < fgkTotUM; i++)
 997     {
 998       for (j = 0; j < fgkRow; j++)
 999         {
1000           cpvTrack[i][j] = new Int_t *[fgkCol];
1001           cpvEdep[i][j]  = new Float_t *[fgkCol];
1002         }
1003     }
1004   for (i = 0; i < fgkTotUM; i++)
1005     {
1006       for (j = 0; j < fgkRow; j++)
1007         {
1008           for (k = 0; k < fgkCol; k++)
1009             {
1010               Int_t nn = fCPVCounter[i][j][k];
1011               if(nn > 0)
1012                 {
1013                   cpvTrack[i][j][k] = new Int_t[nn];
1014                   cpvEdep[i][j][k] = new Float_t[nn];
1015                 }
1016               else
1017                 {
1018                   nn = 1;
1019                   cpvTrack[i][j][k] = new Int_t[nn];
1020                   cpvEdep[i][j][k] = new Float_t[nn];
1021                 }
1022               fCPVCounter[i][j][k] = 0;
1023             }
1024         }
1025     }
1026
1027
1028   Int_t nentries = fCPVCell.GetEntries();
1029
1030   Int_t   mtrackno, ism, ixp, iyp;
1031   Float_t edep;
1032   for (i = 0; i < nentries; i++)
1033     {
1034       AliPMDcell* cpvcell = (AliPMDcell*)fCPVCell.UncheckedAt(i);
1035
1036       mtrackno = cpvcell->GetTrackNumber();
1037       ism      = cpvcell->GetSMNumber();
1038       ixp      = cpvcell->GetX();
1039       iyp      = cpvcell->GetY();
1040       edep     = cpvcell->GetEdep();
1041       Int_t nn = fCPVCounter[ism][ixp][iyp];
1042       cpvTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1043       cpvEdep[ism][ixp][iyp][nn] = edep;
1044       fCPVCounter[ism][ixp][iyp]++;
1045     }
1046
1047   Int_t iz, il;
1048   Int_t im, ix, iy;
1049   Int_t nn;
1050   for (im=0; im<fgkTotUM; im++)
1051     {
1052       for (ix=0; ix<fgkRow; ix++)
1053         {
1054           for (iy=0; iy<fgkCol; iy++)
1055             {
1056               nn = fCPVCounter[im][ix][iy];
1057               if (nn > 1)
1058                 {
1059                   // This block handles if a cell is fired
1060                   // many times by many tracks
1061                   status1  = new Int_t[nn];
1062                   status2  = new Int_t[nn];
1063                   trnarray = new Int_t[nn];
1064                   for (iz = 0; iz < nn; iz++)
1065                     {
1066                       status1[iz] = cpvTrack[im][ix][iy][iz];
1067                     }
1068                   TMath::Sort(nn,status1,status2,jsort);
1069                   Int_t trackOld = -99999;
1070                   Int_t track, trCount = 0;
1071                   for (iz = 0; iz < nn; iz++)
1072                     {
1073                       track = status1[status2[iz]];
1074                       if (trackOld != track)
1075                         {
1076                           trnarray[trCount] = track;
1077                           trCount++;
1078                         }
1079                       trackOld = track;
1080                     }
1081                   delete [] status1;
1082                   delete [] status2;
1083                   Float_t totEdp = 0.;
1084                   trEdp = new Float_t[trCount];
1085                   fracEdp = new Float_t[trCount];
1086                   for (il = 0; il < trCount; il++)
1087                     {
1088                       trEdp[il] = 0.;
1089                       track = trnarray[il];
1090                       for (iz = 0; iz < nn; iz++)
1091                         {
1092                           if (track == cpvTrack[im][ix][iy][iz])
1093                             {
1094                               trEdp[il] += cpvEdep[im][ix][iy][iz];
1095                             }
1096                         }
1097                       totEdp += trEdp[il];
1098                     }
1099                   Int_t ilOld = 0;
1100                   Float_t fracOld = 0.;
1101
1102                   for (il = 0; il < trCount; il++)
1103                     {
1104                       fracEdp[il] = trEdp[il]/totEdp;
1105                       if (fracOld < fracEdp[il])
1106                         {
1107                           fracOld = fracEdp[il];
1108                           ilOld = il;
1109                         }
1110                     }
1111                   fCPVTrackNo[im][ix][iy] = trnarray[ilOld];
1112                   delete [] fracEdp;
1113                   delete [] trEdp;
1114                   delete [] trnarray;
1115                 }
1116               else if (nn == 1)
1117                 {
1118                   // This only handles if a cell is fired
1119                   // by only one track
1120
1121                   fCPVTrackNo[im][ix][iy] = cpvTrack[im][ix][iy][0];
1122
1123                 }
1124               else if (nn ==0)
1125                 {
1126                   // This is if no cell is fired
1127                   fCPVTrackNo[im][ix][iy] = -999;
1128                 }
1129             } // end of iy
1130         } // end of ix
1131     } // end of im
1132
1133   // Delete all the pointers
1134
1135  for (i = 0; i < fgkTotUM; i++)
1136     {
1137       for (j = 0; j < fgkRow; j++)
1138         {
1139           for (k = 0; k < fgkCol; k++)
1140             {
1141               delete []cpvTrack[i][j][k];
1142               delete []cpvEdep[i][j][k];
1143             }
1144         }
1145     }
1146
1147   for (i = 0; i < fgkTotUM; i++)
1148     {
1149       for (j = 0; j < fgkRow; j++)
1150         {
1151           delete [] cpvTrack[i][j];
1152           delete [] cpvEdep[i][j];
1153         }
1154     }
1155
1156   for (i = 0; i < fgkTotUM; i++)
1157     {
1158       delete [] cpvTrack[i];
1159       delete [] cpvEdep[i];
1160     }
1161   delete [] cpvTrack;
1162   delete [] cpvEdep;
1163
1164   //
1165   // End of the cell id assignment
1166   //
1167 }
1168 //____________________________________________________________________________
1169
1170 void AliPMDDigitizer::MergeSDigits(Int_t filenumber, Int_t troffset)
1171 {
1172   // merging sdigits
1173   fRunLoader = AliRunLoader::GetRunLoader(fManager->GetInputFolderName(filenumber));
1174   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
1175   fPMDLoader->LoadSDigits("read");
1176   TTree* treeS = fPMDLoader->TreeS();
1177   AliPMDsdigit  *pmdsdigit;
1178   TBranch *branch = treeS->GetBranch("PMDSDigit");
1179   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1180   branch->SetAddress(&fSDigits);
1181
1182   Int_t   itrackno, itrackpid, idet, ism;
1183   Int_t   ixp, iyp;
1184   Float_t edep;
1185   Int_t nmodules = (Int_t) treeS->GetEntries();
1186   AliDebug(1,Form("Number of Modules in the treeS = %d",nmodules));
1187   AliDebug(1,Form("Track Offset = %d",troffset));
1188   for (Int_t imodule = 0; imodule < nmodules; imodule++)
1189     {
1190       treeS->GetEntry(imodule);
1191       Int_t nentries = fSDigits->GetLast();
1192       AliDebug(2,Form("Number of Entries per Module = %d",nentries));
1193       for (Int_t ient = 0; ient < nentries+1; ient++)
1194         {
1195           pmdsdigit = (AliPMDsdigit*)fSDigits->UncheckedAt(ient);
1196           itrackno  = pmdsdigit->GetTrackNumber();
1197           itrackpid = pmdsdigit->GetTrackPid();
1198           idet      = pmdsdigit->GetDetector();
1199           ism       = pmdsdigit->GetSMNumber();
1200           ixp       = pmdsdigit->GetRow();
1201           iyp       = pmdsdigit->GetColumn();
1202           edep      = pmdsdigit->GetCellEdep();
1203           if (idet == 0)
1204             {
1205               if (fPRE[ism][ixp][iyp] < edep)
1206                 {
1207                   fPRETrackNo[ism][ixp][iyp] = troffset + itrackno;
1208                   fPRETrackPid[ism][ixp][iyp] = itrackpid;
1209                 }
1210               fPRE[ism][ixp][iyp] += edep;
1211             }
1212           else if (idet == 1)
1213             {
1214               if (fCPV[ism][ixp][iyp] < edep)
1215                 {
1216                   fCPVTrackNo[ism][ixp][iyp] = troffset + itrackno;
1217                   fCPVTrackPid[ism][ixp][iyp] = itrackpid;
1218                 }
1219               fCPV[ism][ixp][iyp] += edep;
1220             }
1221         }
1222     }
1223
1224 }
1225 // ----------------------------------------------------------------------
1226 void AliPMDDigitizer::TrackAssignment2Cell()
1227 {
1228   //
1229   // This block assigns the cell id when there are
1230   // multiple tracks in a cell according to the
1231   // energy deposition
1232   //
1233   Bool_t jsort = false;
1234
1235   Int_t i, j, k;
1236
1237   Int_t   *status1;
1238   Int_t   *status2;
1239   Int_t   *trnarray;
1240   Float_t *fracEdp;
1241   Float_t *trEdp;
1242
1243   Int_t   ****pmdTrack;
1244   Float_t ****pmdEdep;
1245
1246   pmdTrack = new Int_t ***[fgkTotUM];
1247   pmdEdep  = new Float_t ***[fgkTotUM];
1248   for (i=0; i<fgkTotUM; i++)
1249     {
1250       pmdTrack[i] = new Int_t **[fgkRow];
1251       pmdEdep[i]  = new Float_t **[fgkRow];
1252     }
1253
1254   for (i = 0; i < fgkTotUM; i++)
1255     {
1256       for (j = 0; j < fgkRow; j++)
1257         {
1258           pmdTrack[i][j] = new Int_t *[fgkCol];
1259           pmdEdep[i][j]  = new Float_t *[fgkCol];
1260         }
1261     }
1262
1263   for (i = 0; i < fgkTotUM; i++)
1264     {
1265       for (j = 0; j < fgkRow; j++)
1266         {
1267           for (k = 0; k < fgkCol; k++)
1268             {
1269               Int_t nn = fPRECounter[i][j][k];
1270               if(nn > 0)
1271                 {
1272                   pmdTrack[i][j][k] = new Int_t[nn];
1273                   pmdEdep[i][j][k] = new Float_t[nn];
1274                 }
1275               else
1276                 {
1277                   nn = 1;
1278                   pmdTrack[i][j][k] = new Int_t[nn];
1279                   pmdEdep[i][j][k] = new Float_t[nn];
1280                 }
1281               fPRECounter[i][j][k] = 0;
1282             }
1283         }
1284     }
1285
1286
1287   Int_t nentries = fCell.GetEntries();
1288
1289   Int_t   mtrackno, ism, ixp, iyp;
1290   Float_t edep;
1291
1292   for (i = 0; i < nentries; i++)
1293     {
1294       AliPMDcell* cell = (AliPMDcell*)fCell.UncheckedAt(i);
1295
1296       mtrackno = cell->GetTrackNumber();
1297       ism      = cell->GetSMNumber();
1298       ixp      = cell->GetX();
1299       iyp      = cell->GetY();
1300       edep     = cell->GetEdep();
1301       Int_t nn = fPRECounter[ism][ixp][iyp];
1302       pmdTrack[ism][ixp][iyp][nn] = (Int_t) mtrackno;
1303       pmdEdep[ism][ixp][iyp][nn] = edep;
1304       fPRECounter[ism][ixp][iyp]++;
1305     }
1306
1307   Int_t iz, il;
1308   Int_t im, ix, iy;
1309   Int_t nn;
1310
1311   for (im=0; im<fgkTotUM; im++)
1312     {
1313       for (ix=0; ix<fgkRow; ix++)
1314         {
1315           for (iy=0; iy<fgkCol; iy++)
1316             {
1317               nn = fPRECounter[im][ix][iy];
1318               if (nn > 1)
1319                 {
1320                   // This block handles if a cell is fired
1321                   // many times by many tracks
1322                   status1  = new Int_t[nn];
1323                   status2  = new Int_t[nn];
1324                   trnarray = new Int_t[nn];
1325                   for (iz = 0; iz < nn; iz++)
1326                     {
1327                       status1[iz] = pmdTrack[im][ix][iy][iz];
1328                     }
1329                   TMath::Sort(nn,status1,status2,jsort);
1330                   Int_t trackOld = -99999;
1331                   Int_t track, trCount = 0;
1332                   for (iz = 0; iz < nn; iz++)
1333                     {
1334                       track = status1[status2[iz]];
1335                       if (trackOld != track)
1336                         {
1337                           trnarray[trCount] = track;
1338                           trCount++;
1339                         }
1340                       trackOld = track;
1341                     }
1342                   delete [] status1;
1343                   delete [] status2;
1344                   Float_t totEdp = 0.;
1345                   trEdp = new Float_t[trCount];
1346                   fracEdp = new Float_t[trCount];
1347                   for (il = 0; il < trCount; il++)
1348                     {
1349                       trEdp[il] = 0.;
1350                       track = trnarray[il];
1351                       for (iz = 0; iz < nn; iz++)
1352                         {
1353                           if (track == pmdTrack[im][ix][iy][iz])
1354                             {
1355                               trEdp[il] += pmdEdep[im][ix][iy][iz];
1356                             }
1357                         }
1358                       totEdp += trEdp[il];
1359                     }
1360                   Int_t ilOld = 0;
1361                   Float_t fracOld = 0.;
1362
1363                   for (il = 0; il < trCount; il++)
1364                     {
1365                       fracEdp[il] = trEdp[il]/totEdp;
1366                       if (fracOld < fracEdp[il])
1367                         {
1368                           fracOld = fracEdp[il];
1369                           ilOld = il;
1370                         }
1371                     }
1372                   fPRETrackNo[im][ix][iy] = trnarray[ilOld];
1373                   delete [] fracEdp;
1374                   delete [] trEdp;
1375                   delete [] trnarray;
1376                 }
1377               else if (nn == 1)
1378                 {
1379                   // This only handles if a cell is fired
1380                   // by only one track
1381
1382                   fPRETrackNo[im][ix][iy] = pmdTrack[im][ix][iy][0];
1383
1384                 }
1385               else if (nn ==0)
1386                 {
1387                   // This is if no cell is fired
1388                   fPRETrackNo[im][ix][iy] = -999;
1389                 }
1390             } // end of iy
1391         } // end of ix
1392     } // end of im
1393
1394   // Delete all the pointers
1395
1396   for (i = 0; i < fgkTotUM; i++)
1397     {
1398       for (j = 0; j < fgkRow; j++)
1399         {
1400           for (k = 0; k < fgkCol; k++)
1401             {
1402               delete [] pmdTrack[i][j][k];
1403               delete [] pmdEdep[i][j][k];
1404             }
1405         }
1406     }
1407
1408   for (i = 0; i < fgkTotUM; i++)
1409     {
1410       for (j = 0; j < fgkRow; j++)
1411         {
1412           delete [] pmdTrack[i][j];
1413           delete [] pmdEdep[i][j];
1414         }
1415     }
1416
1417   for (i = 0; i < fgkTotUM; i++)
1418     {
1419       delete [] pmdTrack[i];
1420       delete [] pmdEdep[i];
1421     }
1422   delete [] pmdTrack;
1423   delete [] pmdEdep;
1424   //
1425   // End of the cell id assignment
1426   //
1427 }
1428 //____________________________________________________________________________
1429 void AliPMDDigitizer::MeV2ADC(Float_t mev, Float_t & adc) const
1430 {
1431   // This converts the simulated edep to ADC according to the
1432   // Test Beam Data
1433   //PS Test in May 2009, Voltage @ 1350 V
1434   // KeV - ADC conversion for 12bit ADC
1435   // MPV data used for the fit and taken here
1436
1437   const Float_t kConstant   = -0.1602;
1438   const Float_t kErConstant =  0.9914;
1439   const Float_t kSlope      = 77.47;
1440   const Float_t kErSlope    =  3.16;
1441
1442   Float_t cons   = gRandom->Gaus(kConstant,kErConstant);
1443   Float_t slop   = gRandom->Gaus(kSlope,kErSlope);
1444
1445   Float_t adc12bit = slop*mev*0.001 + cons;
1446
1447   if (adc12bit < 0.) adc12bit = 0.;
1448
1449   if(adc12bit < 1600.0)
1450     {
1451       adc = (Float_t) adc12bit;
1452     }
1453   else if (adc12bit >= 1600.0)
1454     {
1455       adc = 1600.0;
1456     }
1457 }
1458 //____________________________________________________________________________
1459 void AliPMDDigitizer::AddSDigit(Int_t trnumber, Int_t trpid, Int_t det,
1460                                 Int_t smnumber, Int_t irow, Int_t icol,
1461                                 Float_t adc)
1462 {
1463   // Add SDigit
1464   //
1465   if (!fSDigits) fSDigits = new TClonesArray("AliPMDsdigit", 1000);
1466   TClonesArray &lsdigits = *fSDigits;
1467   new(lsdigits[fNsdigit++])  AliPMDsdigit(trnumber,trpid,det,smnumber,irow,icol,adc);
1468 }
1469 //____________________________________________________________________________
1470
1471 void AliPMDDigitizer::AddDigit(Int_t trnumber, Int_t trpid, Int_t det,
1472                                Int_t smnumber, Int_t irow, Int_t icol,
1473                                Float_t adc)
1474 {
1475   // Add Digit
1476   //
1477   if (!fDigits) fDigits = new TClonesArray("AliPMDdigit", 1000);
1478   TClonesArray &ldigits = *fDigits;
1479   new(ldigits[fNdigit++]) AliPMDdigit(trnumber,trpid, det,smnumber,irow,icol,adc);
1480 }
1481 //____________________________________________________________________________
1482
1483 void AliPMDDigitizer::SetZPosition(Float_t zpos)
1484 {
1485   fZPos = zpos;
1486 }
1487 //____________________________________________________________________________
1488 Float_t AliPMDDigitizer::GetZPosition() const
1489 {
1490   return fZPos;
1491 }
1492 //____________________________________________________________________________
1493
1494 void AliPMDDigitizer::ResetCell()
1495 {
1496   // clears the cell array and also the counter
1497   //  for each cell
1498   //
1499   fCPVCell.Delete();
1500   fCell.Delete();
1501   for (Int_t i = 0; i < fgkTotUM; i++)
1502     {
1503       for (Int_t j = 0; j < fgkRow; j++)
1504         {
1505           for (Int_t k = 0; k < fgkCol; k++)
1506             {
1507               fCPVCounter[i][j][k] = 0;
1508               fPRECounter[i][j][k] = 0;
1509             }
1510         }
1511     }
1512 }
1513 //____________________________________________________________________________
1514 void AliPMDDigitizer::ResetSDigit()
1515 {
1516   // Clears SDigits
1517   fNsdigit = 0;
1518   if (fSDigits) fSDigits->Delete();
1519 }
1520 //____________________________________________________________________________
1521 void AliPMDDigitizer::ResetDigit()
1522 {
1523   // Clears Digits
1524   fNdigit = 0;
1525   if (fDigits) fDigits->Delete();
1526 }
1527 //____________________________________________________________________________
1528
1529 void AliPMDDigitizer::ResetCellADC()
1530 {
1531   // Clears individual cells edep and track number
1532   for (Int_t i = 0; i < fgkTotUM; i++)
1533     {
1534       for (Int_t j = 0; j < fgkRow; j++)
1535         {
1536           for (Int_t k = 0; k < fgkCol; k++)
1537             {
1538               fCPV[i][j][k]         = 0.;
1539               fPRE[i][j][k]         = 0.;
1540               fCPVTrackNo[i][j][k]  = 0;
1541               fPRETrackNo[i][j][k]  = 0;
1542               fCPVTrackPid[i][j][k] = -1;
1543               fPRETrackPid[i][j][k] = -1;
1544             }
1545         }
1546     }
1547 }
1548 //____________________________________________________________________________
1549
1550 void AliPMDDigitizer::UnLoad(Option_t *option)
1551 {
1552   // Unloads all the root files
1553   //
1554   const char *cS = strstr(option,"S");
1555   const char *cD = strstr(option,"D");
1556
1557   fRunLoader->UnloadgAlice();
1558   fRunLoader->UnloadHeader();
1559   fRunLoader->UnloadKinematics();
1560
1561   if (cS)
1562     {
1563       fPMDLoader->UnloadHits();
1564     }
1565   if (cD)
1566     {
1567       fPMDLoader->UnloadHits();
1568       fPMDLoader->UnloadSDigits();
1569     }
1570 }
1571
1572 //----------------------------------------------------------------------
1573 Float_t AliPMDDigitizer::Gain(Int_t det, Int_t smn, Int_t row, Int_t col) const
1574 {
1575   // returns of the gain of the cell
1576   // Added this method by ZA
1577
1578   //cout<<" I am here in gain "<<fCalibData<< "smn,row, col "<<smn
1579   //<<" "<<row<<" "<<col<<endl;
1580
1581   if(!fCalibGain) {
1582     AliError("No calibration data loaded from CDB!!!");
1583     return 1;
1584   }
1585
1586   Float_t GainFact;
1587   GainFact = fCalibGain->GetGainFact(det,smn,row,col);
1588   return GainFact;
1589 }
1590 //----------------------------------------------------------------------
1591 AliPMDCalibData* AliPMDDigitizer::GetCalibGain() const
1592 {
1593   // The run number will be centralized in AliCDBManager,
1594   // you don't need to set it here!
1595   // Added this method by ZA
1596   // Cleaned up by Alberto
1597   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Gain");
1598
1599   if(!entry) AliFatal("Calibration object retrieval failed!");
1600
1601   AliPMDCalibData *calibdata=0;
1602   if (entry) calibdata = (AliPMDCalibData*) entry->GetObject();
1603
1604   if (!calibdata)  AliFatal("No calibration data from calibration database !");
1605
1606   return calibdata;
1607 }
1608 //----------------------------------------------------------------------
1609 AliPMDPedestal* AliPMDDigitizer::GetCalibPed() const
1610 {
1611   // The run number will be centralized in AliCDBManager,
1612   // you don't need to set it here!
1613
1614   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Ped");
1615
1616   if(!entry) AliFatal("Pedestal object retrieval failed!");
1617
1618   AliPMDPedestal *pedestal=0;
1619   if (entry) pedestal = (AliPMDPedestal*) entry->GetObject();
1620
1621   if (!pedestal)  AliFatal("No pedestal data from calibration database !");
1622
1623   return pedestal;
1624 }