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