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