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