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