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