PMD/AliPMDDigitizer.cxx

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