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