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