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