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