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