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