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