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