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