PMD/AliPMDDigitizer.cxx

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