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