PMD/AliPMDClusterFinder.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 //                                                     //
  18 //           Date   : August 05 2003                   //
  19 //  This reads the file PMD.digits.root(TreeD),        //
  20 //  calls the Clustering algorithm and stores the      //
  21 //  clustering output in PMD.RecPoints.root(TreeR)     //
  22 //                                                     //
  23 //-----------------------------------------------------//
  24
  25 #include <Riostream.h>
  26 #include <TTree.h>
  27 #include <TObjArray.h>
  28 #include <TClonesArray.h>
  29
  30 #include "AliLog.h"
  31 #include "AliRunLoader.h"
  32 #include "AliLoader.h"
  33 #include "AliRawReader.h"
  34
  35 #include "AliPMDdigit.h"
  36 #include "AliPMDClusterFinder.h"
  37 #include "AliPMDClustering.h"
  38 #include "AliPMDClusteringV1.h"
  39 #include "AliPMDcluster.h"
  40 #include "AliPMDrecpoint1.h"
  41 #include "AliPMDrechit.h"
  42 #include "AliPMDisocell.h"
  43 #include "AliPMDRawStream.h"
  44 #include "AliPMDCalibData.h"
  45 #include "AliPMDPedestal.h"
  46 #include "AliPMDddldata.h"
  47
  48 #include "AliDAQ.h"
  49 #include "AliCDBManager.h"
  50 #include "AliCDBEntry.h"
  51
  52
  53
  54 ClassImp(AliPMDClusterFinder)
  55
  56 AliPMDClusterFinder::AliPMDClusterFinder():
  57   fRunLoader(0),
  58   fPMDLoader(0),
  59   fCalibGain(GetCalibGain()),
  60   fCalibPed(GetCalibPed()),
  61   fTreeD(0),
  62   fTreeR(0),
  63   fDigits(new TClonesArray("AliPMDdigit", 1000)),
  64   fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
  65   fRechits(new TClonesArray("AliPMDrechit", 1000)),
  66   fNpoint(0),
  67   fNhit(0),
  68   fDetNo(0),
  69   fEcut(0.)
  70 {
  71 //
  72 // Constructor
  73 //
  74 }
  75 // ------------------------------------------------------------------------- //
  76 AliPMDClusterFinder::AliPMDClusterFinder(AliRunLoader* runLoader):
  77   fRunLoader(runLoader),
  78   fPMDLoader(runLoader->GetLoader("PMDLoader")),
  79   fCalibGain(GetCalibGain()),
  80   fCalibPed(GetCalibPed()),
  81   fTreeD(0),
  82   fTreeR(0),
  83   fDigits(new TClonesArray("AliPMDdigit", 1000)),
  84   fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
  85   fRechits(new TClonesArray("AliPMDrechit", 1000)),
  86   fNpoint(0),
  87   fNhit(0),
  88   fDetNo(0),
  89   fEcut(0.)
  90 {
  91 //
  92 // Constructor
  93 //
  94 }
  95 // ------------------------------------------------------------------------- //
  96 AliPMDClusterFinder::AliPMDClusterFinder(const AliPMDClusterFinder & finder):
  97   TObject(finder),
  98   fRunLoader(0),
  99   fPMDLoader(0),
 100   fCalibGain(GetCalibGain()),
 101   fCalibPed(GetCalibPed()),
 102   fTreeD(0),
 103   fTreeR(0),
 104   fDigits(NULL),
 105   fRecpoints(NULL),
 106   fRechits(NULL),
 107   fNpoint(0),
 108   fNhit(0),
 109   fDetNo(0),
 110   fEcut(0.)
 111 {
 112   // copy constructor
 113   AliError("Copy constructor not allowed");
 114 }
 115 // ------------------------------------------------------------------------- //
 116 AliPMDClusterFinder &AliPMDClusterFinder::operator=(const AliPMDClusterFinder & /*finder*/)
 117 {
 118  // assignment op
 119   AliError("Assignment Operator not allowed");
 120   return *this;
 121 }
 122 // ------------------------------------------------------------------------- //
 123 AliPMDClusterFinder::~AliPMDClusterFinder()
 124 {
 125   // Destructor
 126   if (fDigits)
 127     {
 128         fDigits->Clear();
 129     }
 130   if (fRecpoints)
 131     {
 132       fRecpoints->Clear();
 133     }
 134   if (fRechits)
 135     {
 136       fRechits->Clear();
 137     }
 138
 139 }
 140 // ------------------------------------------------------------------------- //
 141
 142 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt)
 143 {
 144   // Converts digits to recpoints after running clustering
 145   // algorithm on CPV plane and PREshower plane
 146   //
 147
 148   Int_t    det = 0,smn = 0;
 149   Int_t    xpos,ypos;
 150   Float_t  adc;
 151   Int_t    ismn;
 152   Int_t    idet;
 153   Float_t  clusdata[6];
 154
 155   AliPMDisocell *pmdiso = 0x0;
 156
 157   TObjArray *pmdcont = new TObjArray();
 158   TObjArray *pmdisocell = new TObjArray();
 159
 160   AliPMDClustering *pmdclust = new AliPMDClusteringV1();
 161
 162   pmdclust->SetEdepCut(fEcut);
 163
 164   fRunLoader->GetEvent(ievt);
 165
 166
 167   fTreeD = fPMDLoader->TreeD();
 168   if (fTreeD == 0x0)
 169     {
 170       AliFatal("AliPMDClusterFinder: Can not get TreeD");
 171
 172     }
 173   AliPMDdigit  *pmddigit;
 174   TBranch *branch = fTreeD->GetBranch("PMDDigit");
 175   branch->SetAddress(&fDigits);
 176
 177   ResetRecpoint();
 178
 179   fTreeR = fPMDLoader->TreeR();
 180   if (fTreeR == 0x0)
 181     {
 182       fPMDLoader->MakeTree("R");
 183       fTreeR = fPMDLoader->TreeR();
 184     }
 185
 186   Int_t bufsize = 16000;
 187   TBranch * branch1 = fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize);
 188   TBranch * branch2 = fTreeR->Branch("PMDRechit", &fRechits, bufsize);
 189
 190   Int_t nmodules = (Int_t) fTreeD->GetEntries();
 191
 192   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 193     {
 194       ResetCellADC();
 195       fTreeD->GetEntry(imodule);
 196       Int_t nentries = fDigits->GetLast();
 197       for (Int_t ient = 0; ient < nentries+1; ient++)
 198         {
 199           pmddigit = (AliPMDdigit*)fDigits->UncheckedAt(ient);
 200
 201           det    = pmddigit->GetDetector();
 202           smn    = pmddigit->GetSMNumber();
 203           xpos   = pmddigit->GetRow();
 204           ypos   = pmddigit->GetColumn();
 205           adc    = pmddigit->GetADC();
 206           if(xpos < 0 || xpos > 48 || ypos < 0 || ypos > 96)
 207             {
 208               AliError(Form("*Row %d and Column NUMBER %d NOT Valid *",
 209                               xpos, ypos));
 210               continue;
 211             }
 212           // CALIBRATION
 213           Float_t gain = fCalibGain->GetGainFact(det,smn,xpos,ypos);
 214           // printf("adc = %d gain = %f\n",adc,gain);
 215
 216           adc = adc*gain;
 217
 218           //Int_t trno   = pmddigit->GetTrackNumber();
 219           fCellTrack[xpos][ypos] = pmddigit->GetTrackNumber();
 220           fCellPid[xpos][ypos]   = pmddigit->GetTrackPid();
 221           fCellADC[xpos][ypos]   = (Double_t) adc;
 222         }
 223
 224       idet = det;
 225       ismn = smn;
 226       //pmdclust->DoClust(idet,ismn,fCellADC,pmdisocell,pmdcont);
 227       pmdclust->DoClust(idet,ismn,fCellTrack,fCellPid,fCellADC,
 228                         pmdisocell,pmdcont);
 229
 230       Int_t nentries1 = pmdcont->GetEntries();
 231
 232       AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 233
 234       for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 235         {
 236           AliPMDcluster *pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
 237           idet        = pmdcl->GetDetector();
 238           ismn        = pmdcl->GetSMN();
 239           clusdata[0] = pmdcl->GetClusX();
 240           clusdata[1] = pmdcl->GetClusY();
 241           clusdata[2] = pmdcl->GetClusADC();
 242           clusdata[3] = pmdcl->GetClusCells();
 243           clusdata[4] = pmdcl->GetClusSigmaX();
 244           clusdata[5] = pmdcl->GetClusSigmaY();
 245
 246           AddRecPoint(idet,ismn,clusdata);
 247
 248           Int_t ncell = (Int_t) clusdata[3];
 249           for(Int_t ihit = 0; ihit < ncell; ihit++)
 250             {
 251               Int_t celldataX = pmdcl->GetClusCellX(ihit);
 252               Int_t celldataY = pmdcl->GetClusCellY(ihit);
 253               Int_t celldataTr = pmdcl->GetClusCellTrack(ihit);
 254               Int_t celldataPid   = pmdcl->GetClusCellPid(ihit);
 255               AddRecHit(celldataX, celldataY, celldataTr, celldataPid);
 256               //AddRecHit(celldataX, celldataY);
 257             }
 258           branch2->Fill();
 259           ResetRechit();
 260         }
 261       pmdcont->Delete();
 262
 263           // Added single isolated cell for offline gain calibration
 264           nentries1 = pmdisocell->GetEntries();
 265           AliDebug(1,Form("Total number of isolated single cell clusters = %d",nentries1));
 266
 267           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 268             {
 269               pmdiso = (AliPMDisocell*)pmdisocell->UncheckedAt(ient1);
 270               idet = pmdiso->GetDetector();
 271               ismn = pmdiso->GetSmn();
 272               clusdata[0] = (Float_t) pmdiso->GetRow();
 273               clusdata[1] = (Float_t) pmdiso->GetCol();
 274               clusdata[2] = pmdiso->GetADC();
 275               clusdata[3] = 1.;
 276               clusdata[4] = -99.;
 277               clusdata[5] = -99.;
 278
 279               AddRecPoint(idet,ismn,clusdata);
 280             }
 281           pmdisocell->Delete();
 282
 283       branch1->Fill();
 284       ResetRecpoint();
 285
 286     } // modules
 287
 288   ResetCellADC();
 289   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 290   fPMDLoader->WriteRecPoints("OVERWRITE");
 291
 292   //   delete the pointers
 293   delete pmdclust;
 294   delete pmdcont;
 295   delete pmdisocell;
 296
 297 }
 298 // ------------------------------------------------------------------------- //
 299
 300 void AliPMDClusterFinder::Digits2RecPoints(TTree *digitsTree,
 301                                            TTree *clustersTree)
 302 {
 303   // Converts digits to recpoints after running clustering
 304   // algorithm on CPV plane and PREshower plane
 305   //
 306   // This algorithm is called during the reconstruction from digits
 307
 308   Int_t    det = 0,smn = 0;
 309   Int_t    xpos,ypos;
 310   Float_t  adc;
 311   Int_t    ismn;
 312   Int_t    idet;
 313   Float_t  clusdata[6];
 314
 315   AliPMDcluster *pmdcl = 0x0;
 316   AliPMDisocell *pmdiso = 0x0;
 317
 318   TObjArray *pmdcont = new TObjArray();
 319   TObjArray *pmdisocell = new TObjArray();
 320   AliPMDClustering *pmdclust = new AliPMDClusteringV1();
 321
 322   pmdclust->SetEdepCut(fEcut);
 323
 324   AliPMDdigit  *pmddigit;
 325   TBranch *branch = digitsTree->GetBranch("PMDDigit");
 326   branch->SetAddress(&fDigits);
 327
 328   ResetRecpoint();
 329
 330   Int_t bufsize = 16000;
 331   TBranch * branch1 = clustersTree->Branch("PMDRecpoint", &fRecpoints, bufsize);
 332   TBranch * branch2 = clustersTree->Branch("PMDRechit", &fRechits, bufsize);
 333
 334   Int_t nmodules = (Int_t) digitsTree->GetEntries();
 335
 336   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 337     {
 338
 339       Int_t totADCMod = 0;
 340       ResetCellADC();
 341       digitsTree->GetEntry(imodule);
 342       Int_t nentries = fDigits->GetLast();
 343       for (Int_t ient = 0; ient < nentries+1; ient++)
 344         {
 345           pmddigit = (AliPMDdigit*)fDigits->UncheckedAt(ient);
 346
 347           det    = pmddigit->GetDetector();
 348           smn    = pmddigit->GetSMNumber();
 349           xpos   = pmddigit->GetRow();
 350           ypos   = pmddigit->GetColumn();
 351           adc    = pmddigit->GetADC();
 352           if(xpos < 0 || xpos > 48 || ypos < 0 || ypos > 96)
 353             {
 354               AliError(Form("*Row %d and Column NUMBER %d NOT Valid *",
 355                             xpos, ypos));
 356               continue;
 357             }
 358
 359           // Pedestal Subtraction
 360           Int_t   pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,xpos,ypos);
 361           Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 362           Float_t pedrms     = (Float_t)pedrms1/10.;
 363           Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 364           //printf("%f %f\n",pedmean, pedrms);
 365
 366           Float_t adc1 = adc - (pedmean + 3.0*pedrms);
 367
 368           // CALIBRATION
 369           Float_t gain = fCalibGain->GetGainFact(det,smn,xpos,ypos);
 370           // printf("adc = %d gain = %f\n",adc,gain);
 371
 372           adc = adc1*gain;
 373
 374           //Int_t trno   = pmddigit->GetTrackNumber();
 375           fCellTrack[xpos][ypos] = pmddigit->GetTrackNumber();
 376           fCellPid[xpos][ypos] = pmddigit->GetTrackPid();
 377           fCellADC[xpos][ypos] = (Double_t) adc;
 378
 379           totADCMod += (Int_t) adc;
 380
 381         }
 382
 383       idet = det;
 384       ismn = smn;
 385
 386       if (totADCMod <= 0) continue;
 387
 388       pmdclust->DoClust(idet,ismn,fCellTrack,fCellPid,fCellADC,
 389                         pmdisocell,pmdcont);
 390
 391       Int_t nentries1 = pmdcont->GetEntries();
 392
 393       AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 394
 395       for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 396         {
 397           pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
 398           idet        = pmdcl->GetDetector();
 399           ismn        = pmdcl->GetSMN();
 400           clusdata[0] = pmdcl->GetClusX();
 401           clusdata[1] = pmdcl->GetClusY();
 402           clusdata[2] = pmdcl->GetClusADC();
 403           clusdata[3] = pmdcl->GetClusCells();
 404           clusdata[4] = pmdcl->GetClusSigmaX();
 405           clusdata[5] = pmdcl->GetClusSigmaY();
 406
 407           AddRecPoint(idet,ismn,clusdata);
 408
 409           Int_t ncell = (Int_t) clusdata[3];
 410           for(Int_t ihit = 0; ihit < ncell; ihit++)
 411             {
 412               Int_t celldataX = pmdcl->GetClusCellX(ihit);
 413               Int_t celldataY = pmdcl->GetClusCellY(ihit);
 414               Int_t celldataTr = pmdcl->GetClusCellTrack(ihit);
 415               Int_t celldataPid   = pmdcl->GetClusCellPid(ihit);
 416               AddRecHit(celldataX, celldataY, celldataTr, celldataPid);
 417             }
 418           branch2->Fill();
 419           ResetRechit();
 420         }
 421       pmdcont->Delete();
 422
 423       // Added single isolated cell for offline gain calibration
 424       nentries1 = pmdisocell->GetEntries();
 425       AliDebug(1,Form("Total number of isolated single cell clusters = %d",nentries1));
 426
 427       for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 428         {
 429           pmdiso = (AliPMDisocell*)pmdisocell->UncheckedAt(ient1);
 430           idet = pmdiso->GetDetector();
 431           ismn = pmdiso->GetSmn();
 432           clusdata[0] = (Float_t) pmdiso->GetRow();
 433           clusdata[1] = (Float_t) pmdiso->GetCol();
 434           clusdata[2] = pmdiso->GetADC();
 435           clusdata[3] = 1.;
 436           clusdata[4] = -99.;
 437           clusdata[5] = -99.;
 438
 439           AddRecPoint(idet,ismn,clusdata);
 440         }
 441       pmdisocell->Delete();
 442
 443       branch1->Fill();
 444       ResetRecpoint();
 445
 446     } // modules
 447
 448
 449   ResetCellADC();
 450
 451   //   delete the pointers
 452   delete pmdclust;
 453   delete pmdcont;
 454   delete pmdisocell;
 455 }
 456 // ------------------------------------------------------------------------- //
 457
 458 void AliPMDClusterFinder::Digits2RecPoints(AliRawReader *rawReader,
 459                                            TTree *clustersTree)
 460 {
 461   // Converts RAW data to recpoints after running clustering
 462   // algorithm on CPV and PREshower plane
 463   //
 464   // This method is called at the time of reconstruction from RAW data
 465
 466
 467   AliPMDddldata *pmdddl = 0x0;
 468   AliPMDcluster *pmdcl  = 0x0;
 469   AliPMDisocell *pmdiso = 0x0;
 470
 471
 472   Float_t  clusdata[6];
 473   TObjArray pmdddlcont;
 474
 475   TObjArray *pmdcont = new TObjArray();
 476   TObjArray *pmdisocell = new TObjArray();
 477
 478   AliPMDClustering *pmdclust = new AliPMDClusteringV1();
 479
 480   pmdclust->SetEdepCut(fEcut);
 481
 482   ResetRecpoint();
 483
 484   Int_t bufsize = 16000;
 485   TBranch *branch1 = clustersTree->Branch("PMDRecpoint", &fRecpoints, bufsize);
 486
 487   TBranch * branch2 = clustersTree->Branch("PMDRechit", &fRechits, bufsize);
 488
 489   const Int_t kRow = 48;
 490   const Int_t kCol = 96;
 491
 492   Int_t idet = 0;
 493   Int_t iSMN = 0;
 494
 495   Int_t indexDDL = -1;
 496   AliPMDRawStream pmdinput(rawReader);
 497
 498   while ((indexDDL = pmdinput.DdlData(&pmdddlcont)) >=0)
 499   {
 500       if (indexDDL < 4)
 501         {
 502           iSMN = 6;
 503         }
 504       else if (indexDDL >= 4)
 505         {
 506           iSMN = 12;
 507         }
 508       Int_t ***precpvADC;
 509       precpvADC = new int **[iSMN];
 510       for (Int_t i=0; i<iSMN; i++) precpvADC[i] = new int *[kRow];
 511       for (Int_t i=0; i<iSMN;i++)
 512         {
 513           for (Int_t j=0; j<kRow; j++) precpvADC[i][j] = new int [kCol];
 514         }
 515       for (Int_t i = 0; i < iSMN; i++)
 516         {
 517           for (Int_t j = 0; j < kRow; j++)
 518             {
 519               for (Int_t k = 0; k < kCol; k++)
 520                 {
 521                   precpvADC[i][j][k] = 0;
 522                 }
 523             }
 524         }
 525       ResetCellADC();
 526
 527       Int_t indexsmn = 0;
 528       Int_t ientries = pmdddlcont.GetEntries();
 529       for (Int_t ient = 0; ient < ientries; ient++)
 530         {
 531           pmdddl = (AliPMDddldata*)pmdddlcont.UncheckedAt(ient);
 532
 533           Int_t det = pmdddl->GetDetector();
 534           Int_t smn = pmdddl->GetSMN();
 535           //Int_t mcm = pmdddl->GetMCM();
 536           //Int_t chno = pmdddl->GetChannel();
 537           Int_t row = pmdddl->GetRow();
 538           Int_t col = pmdddl->GetColumn();
 539           Int_t sig = pmdddl->GetSignal();
 540
 541           if(smn == -1)
 542             {
 543               AliError(Form("*MODULE NUMBER WRONG %d *",smn));
 544               continue;
 545             }
 546           if(row < 0 || row > 48 || col < 0 || col > 96)
 547             {
 548               AliError(Form("*Row %d and Column NUMBER %d NOT Valid *",
 549                             row, col));
 550               continue;
 551             }
 552
 553           // Pedestal Subtraction
 554           Int_t   pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,row,col);
 555           Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 556           Float_t pedrms     = (Float_t)pedrms1/10.;
 557           Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 558
 559           //printf("%f %f\n",pedmean, pedrms);
 560
 561           // Float_t sig1 = (Float_t) sig;
 562           Float_t sig1 = (Float_t) sig - (pedmean + 3.0*pedrms);
 563
 564           // CALIBRATION
 565           Float_t gain = fCalibGain->GetGainFact(det,smn,row,col);
 566           //printf("sig = %d gain = %f\n",sig,gain);
 567           sig = (Int_t) (sig1*gain);
 568
 569           if (indexDDL < 4)
 570             {
 571               if (det != 0)
 572                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 573                               indexDDL, det));
 574               indexsmn = smn - indexDDL * 6;
 575             }
 576           else if (indexDDL == 4)
 577             {
 578               if (det != 1)
 579                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 580                               indexDDL, det));
 581               if (smn < 6)
 582                 {
 583                   indexsmn = smn;
 584                 }
 585               else if (smn >= 18 && smn < 24)
 586                 {
 587                   indexsmn = smn - 12;
 588                 }
 589             }
 590           else if (indexDDL == 5)
 591             {
 592               if (det != 1)
 593                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 594                               indexDDL, det));
 595               if (smn >= 6 && smn < 12)
 596                 {
 597                   indexsmn = smn - 6;
 598                 }
 599               else if (smn >= 12 && smn < 18)
 600                 {
 601                   indexsmn = smn - 6;
 602                 }
 603             }
 604           precpvADC[indexsmn][row][col] = sig;
 605         }
 606
 607       pmdddlcont.Delete();
 608
 609       Int_t totAdcMod = 0;
 610
 611       Int_t ismn = 0;
 612       for (indexsmn = 0; indexsmn < iSMN; indexsmn++)
 613         {
 614           ResetCellADC();
 615           totAdcMod = 0;
 616           for (Int_t irow = 0; irow < kRow; irow++)
 617             {
 618               for (Int_t icol = 0; icol < kCol; icol++)
 619                 {
 620                   fCellTrack[irow][icol] = -1;
 621                   fCellPid[irow][icol]   = -1;
 622
 623                   fCellADC[irow][icol] =
 624                     (Double_t) precpvADC[indexsmn][irow][icol];
 625                   totAdcMod += precpvADC[indexsmn][irow][icol];
 626                 } // row
 627             }     // col
 628
 629           if (indexDDL < 4)
 630             {
 631               ismn = indexsmn + indexDDL * 6;
 632               idet = 0;
 633             }
 634           else if (indexDDL == 4)
 635             {
 636               if (indexsmn < 6)
 637                 {
 638                   ismn = indexsmn;
 639                 }
 640               else if (indexsmn >= 6 && indexsmn < 12)
 641                 {
 642                   ismn = indexsmn + 12;
 643                 }
 644               idet = 1;
 645             }
 646           else if (indexDDL == 5)
 647             {
 648               if (indexsmn < 6)
 649                 {
 650                   ismn = indexsmn + 6;
 651                 }
 652               else if (indexsmn >= 6 && indexsmn < 12)
 653                 {
 654                   ismn = indexsmn + 6;
 655                 }
 656               idet = 1;
 657             }
 658
 659           if (totAdcMod <= 0) continue;
 660
 661
 662           //pmdclust->DoClust(idet,ismn,fCellADC,pmdisocell,pmdcont);
 663           pmdclust->DoClust(idet,ismn,fCellTrack,fCellPid,fCellADC,
 664                             pmdisocell,pmdcont);
 665
 666           Int_t nentries1 = pmdcont->GetEntries();
 667
 668           AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 669
 670           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 671             {
 672               pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
 673               idet        = pmdcl->GetDetector();
 674               ismn        = pmdcl->GetSMN();
 675               clusdata[0] = pmdcl->GetClusX();
 676               clusdata[1] = pmdcl->GetClusY();
 677               clusdata[2] = pmdcl->GetClusADC();
 678               clusdata[3] = pmdcl->GetClusCells();
 679               clusdata[4] = pmdcl->GetClusSigmaX();
 680               clusdata[5] = pmdcl->GetClusSigmaY();
 681
 682               AddRecPoint(idet,ismn,clusdata);
 683
 684               Int_t ncell = (Int_t) clusdata[3];
 685               for(Int_t ihit = 0; ihit < ncell; ihit++)
 686                 {
 687                   Int_t celldataX = pmdcl->GetClusCellX(ihit);
 688                   Int_t celldataY = pmdcl->GetClusCellY(ihit);
 689                   Int_t celldataTr = pmdcl->GetClusCellTrack(ihit);
 690                   Int_t celldataPid   = pmdcl->GetClusCellPid(ihit);
 691                   AddRecHit(celldataX, celldataY, celldataTr, celldataPid);
 692                   //AddRecHit(celldataX, celldataY);
 693                 }
 694               branch2->Fill();
 695               ResetRechit();
 696
 697             }
 698           pmdcont->Delete();
 699
 700
 701           // Added single isolated cell for offline gain calibration
 702           nentries1 = pmdisocell->GetEntries();
 703           AliDebug(1,Form("Total number of isolated single cell clusters = %d",nentries1));
 704
 705           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 706             {
 707               pmdiso = (AliPMDisocell*)pmdisocell->UncheckedAt(ient1);
 708               idet = pmdiso->GetDetector();
 709               ismn = pmdiso->GetSmn();
 710               clusdata[0] = (Float_t) pmdiso->GetRow();
 711               clusdata[1] = (Float_t) pmdiso->GetCol();
 712               clusdata[2] = pmdiso->GetADC();
 713               clusdata[3] = 1.;
 714               clusdata[4] = -99.;
 715               clusdata[5] = -99.;
 716
 717               AddRecPoint(idet,ismn,clusdata);
 718             }
 719           pmdisocell->Delete();
 720
 721           branch1->Fill();
 722           ResetRecpoint();
 723
 724
 725         } // smn
 726
 727       for (Int_t i=0; i<iSMN; i++)
 728         {
 729           for (Int_t j=0; j<kRow; j++) delete [] precpvADC[i][j];
 730         }
 731       for (Int_t i=0; i<iSMN; i++) delete [] precpvADC[i];
 732       delete [] precpvADC;
 733
 734     } // DDL Loop
 735
 736
 737   ResetCellADC();
 738
 739   //   delete the pointers
 740   delete pmdclust;
 741   delete pmdcont;
 742   delete pmdisocell;
 743
 744 }
 745 // ------------------------------------------------------------------------- //
 746
 747 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt, AliRawReader *rawReader)
 748 {
 749   // Converts RAW data to recpoints after running clustering
 750   // algorithm on CPV and PREshower plane
 751   //
 752
 753   Float_t  clusdata[6];
 754
 755   TObjArray pmdddlcont;
 756
 757   //AliPMDcluster *pmdcl  = 0x0;
 758   AliPMDisocell *pmdiso  = 0x0;
 759
 760
 761   TObjArray *pmdcont = new TObjArray();
 762   TObjArray *pmdisocell = new TObjArray();
 763
 764   AliPMDClustering *pmdclust = new AliPMDClusteringV1();
 765
 766   pmdclust->SetEdepCut(fEcut);
 767
 768   fRunLoader->GetEvent(ievt);
 769
 770   ResetRecpoint();
 771
 772   fTreeR = fPMDLoader->TreeR();
 773   if (fTreeR == 0x0)
 774     {
 775       fPMDLoader->MakeTree("R");
 776       fTreeR = fPMDLoader->TreeR();
 777     }
 778   Int_t bufsize = 16000;
 779   TBranch *branch1 = fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize);
 780   TBranch *branch2 = fTreeR->Branch("PMDRechit", &fRechits, bufsize);
 781
 782   const Int_t kRow = 48;
 783   const Int_t kCol = 96;
 784
 785   Int_t idet = 0;
 786   Int_t iSMN = 0;
 787
 788   AliPMDRawStream pmdinput(rawReader);
 789   Int_t indexDDL = -1;
 790
 791   while ((indexDDL = pmdinput.DdlData(&pmdddlcont)) >=0) {
 792
 793       if (indexDDL < 4)
 794         {
 795           iSMN = 6;
 796         }
 797       else if (indexDDL >= 4)
 798         {
 799           iSMN = 12;
 800         }
 801       Int_t ***precpvADC;
 802       precpvADC = new int **[iSMN];
 803       for (Int_t i=0; i<iSMN; i++) precpvADC[i] = new int *[kRow];
 804       for (Int_t i=0; i<iSMN;i++)
 805         {
 806           for (Int_t j=0; j<kRow; j++) precpvADC[i][j] = new int [kCol];
 807         }
 808       for (Int_t i = 0; i < iSMN; i++)
 809         {
 810           for (Int_t j = 0; j < kRow; j++)
 811             {
 812               for (Int_t k = 0; k < kCol; k++)
 813                 {
 814                   precpvADC[i][j][k] = 0;
 815                 }
 816             }
 817         }
 818       ResetCellADC();
 819
 820
 821       Int_t indexsmn = 0;
 822       Int_t ientries = pmdddlcont.GetEntries();
 823       for (Int_t ient = 0; ient < ientries; ient++)
 824         {
 825           AliPMDddldata *pmdddl = (AliPMDddldata*)pmdddlcont.UncheckedAt(ient);
 826
 827           Int_t det = pmdddl->GetDetector();
 828           Int_t smn = pmdddl->GetSMN();
 829           //Int_t mcm = pmdddl->GetMCM();
 830           //Int_t chno = pmdddl->GetChannel();
 831           Int_t row = pmdddl->GetRow();
 832           Int_t col = pmdddl->GetColumn();
 833           Int_t sig = pmdddl->GetSignal();
 834           if(row < 0 || row > 48 || col < 0 || col > 96)
 835             {
 836               AliError(Form("*Row %d and Column NUMBER %d NOT Valid *",
 837                               row, col));
 838               continue;
 839             }
 840           // Pedestal Subtraction
 841           Int_t   pedmeanrms = fCalibPed->GetPedMeanRms(det,smn,row,col);
 842           Int_t   pedrms1    = (Int_t) pedmeanrms%100;
 843           Float_t pedrms     = (Float_t)pedrms1/10.;
 844           Float_t pedmean    = (Float_t) (pedmeanrms - pedrms1)/1000.0;
 845
 846           //printf("%f %f\n",pedmean, pedrms);
 847
 848           //Float_t sig1 = (Float_t) sig;
 849           Float_t sig1 = (Float_t) sig - (pedmean + 3.0*pedrms);
 850           // CALIBRATION
 851           Float_t gain = fCalibGain->GetGainFact(det,smn,row,col);
 852
 853           //printf("sig = %d gain = %f\n",sig,gain);
 854           sig = (Int_t) (sig1*gain);
 855
 856
 857           if (indexDDL < 4)
 858             {
 859               if (det != 0)
 860                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 861                               indexDDL, det));
 862               indexsmn = smn - indexDDL * 6;
 863             }
 864           else if (indexDDL == 4)
 865             {
 866               if (det != 1)
 867                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 868                               indexDDL, det));
 869               if (smn < 6)
 870                 {
 871                   indexsmn = smn;
 872                 }
 873               else if (smn >= 18 && smn < 24)
 874                 {
 875                   indexsmn = smn - 12;
 876                 }
 877             }
 878           else if (indexDDL == 5)
 879             {
 880               if (det != 1)
 881                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 882                               indexDDL, det));
 883               if (smn >= 6 && smn < 12)
 884                 {
 885                   indexsmn = smn - 6;
 886                 }
 887               else if (smn >= 12 && smn < 18)
 888                 {
 889                   indexsmn = smn - 6;
 890                 }
 891             }
 892           precpvADC[indexsmn][row][col] = sig;
 893
 894         }
 895
 896       pmdddlcont.Delete();
 897
 898       Int_t ismn = 0;
 899       for (indexsmn = 0; indexsmn < iSMN; indexsmn++)
 900         {
 901           ResetCellADC();
 902           for (Int_t irow = 0; irow < kRow; irow++)
 903             {
 904               for (Int_t icol = 0; icol < kCol; icol++)
 905                 {
 906                   fCellTrack[irow][icol] = -1;
 907                   fCellPid[irow][icol]   = -1;
 908                   fCellADC[irow][icol] =
 909                     (Double_t) precpvADC[indexsmn][irow][icol];
 910                 } // row
 911             }     // col
 912
 913
 914           if (indexDDL < 4)
 915             {
 916               ismn = indexsmn + indexDDL * 6;
 917               idet = 0;
 918             }
 919           else if (indexDDL == 4)
 920             {
 921               if (indexsmn < 6)
 922                 {
 923                   ismn = indexsmn;
 924                 }
 925               else if (indexsmn >= 6 && indexsmn < 12)
 926                 {
 927                   ismn = indexsmn + 12;
 928                 }
 929               idet = 1;
 930             }
 931           else if (indexDDL == 5)
 932             {
 933               if (indexsmn < 6)
 934                 {
 935                   ismn = indexsmn + 6;
 936                 }
 937               else if (indexsmn >= 6 && indexsmn < 12)
 938                 {
 939                   ismn = indexsmn + 6;
 940                 }
 941               idet = 1;
 942             }
 943
 944           //pmdclust->DoClust(idet,ismn,fCellADC,pmdisocell,pmdcont);
 945           pmdclust->DoClust(idet,ismn,fCellTrack,fCellPid,fCellADC,
 946                             pmdisocell,pmdcont);
 947
 948           Int_t nentries1 = pmdcont->GetEntries();
 949
 950           AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 951
 952           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 953             {
 954               AliPMDcluster *pmdcl =
 955                 (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
 956               idet        = pmdcl->GetDetector();
 957               ismn        = pmdcl->GetSMN();
 958               clusdata[0] = pmdcl->GetClusX();
 959               clusdata[1] = pmdcl->GetClusY();
 960               clusdata[2] = pmdcl->GetClusADC();
 961               clusdata[3] = pmdcl->GetClusCells();
 962               clusdata[4] = pmdcl->GetClusSigmaX();
 963               clusdata[5] = pmdcl->GetClusSigmaY();
 964
 965               AddRecPoint(idet,ismn,clusdata);
 966
 967               Int_t ncell = (Int_t) clusdata[3];
 968               for(Int_t ihit = 0; ihit < ncell; ihit++)
 969                 {
 970                   Int_t celldataX = pmdcl->GetClusCellX(ihit);
 971                   Int_t celldataY = pmdcl->GetClusCellY(ihit);
 972                   Int_t celldataTr = pmdcl->GetClusCellTrack(ihit);
 973                   Int_t celldataPid = pmdcl->GetClusCellPid(ihit);
 974                   AddRecHit(celldataX, celldataY, celldataTr, celldataPid);
 975                   //AddRecHit(celldataX, celldataY);
 976                 }
 977               branch2->Fill();
 978               ResetRechit();
 979
 980             }
 981           pmdcont->Delete();
 982
 983           // Added single isolated cell for offline gain calibration
 984           nentries1 = pmdisocell->GetEntries();
 985           AliDebug(1,Form("Total number of isolated single cell clusters = %d",nentries1));
 986
 987           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 988             {
 989               pmdiso = (AliPMDisocell*)pmdisocell->UncheckedAt(ient1);
 990               idet = pmdiso->GetDetector();
 991               ismn = pmdiso->GetSmn();
 992               clusdata[0] = (Float_t) pmdiso->GetRow();
 993               clusdata[1] = (Float_t) pmdiso->GetCol();
 994               clusdata[2] = pmdiso->GetADC();
 995               clusdata[3] = 1.;
 996               clusdata[4] = -99.;
 997               clusdata[5] = -99.;
 998
 999               AddRecPoint(idet,ismn,clusdata);
1000             }
1001           pmdisocell->Delete();
1002
1003           branch1->Fill();
1004           ResetRecpoint();
1005
1006
1007         } // smn
1008
1009       for (Int_t i=0; i<iSMN; i++)
1010         {
1011           for (Int_t j=0; j<kRow; j++) delete [] precpvADC[i][j];
1012         }
1013       for (Int_t i=0; i<iSMN; i++) delete [] precpvADC[i];
1014       delete precpvADC;
1015     } // DDL Loop
1016
1017
1018   ResetCellADC();
1019
1020   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
1021   fPMDLoader->WriteRecPoints("OVERWRITE");
1022
1023   //   delete the pointers
1024   delete pmdclust;
1025   delete pmdcont;
1026   delete pmdisocell;
1027 }
1028 // ------------------------------------------------------------------------- //
1029 void AliPMDClusterFinder::SetCellEdepCut(Float_t ecut)
1030 {
1031   fEcut = ecut;
1032 }
1033 // ------------------------------------------------------------------------- //
1034 void AliPMDClusterFinder::AddRecPoint(Int_t idet,Int_t ismn,Float_t *clusdata)
1035 {
1036   // Add Reconstructed points
1037   //
1038   TClonesArray &lrecpoints = *fRecpoints;
1039   AliPMDrecpoint1 *newrecpoint;
1040   newrecpoint = new AliPMDrecpoint1(idet, ismn, clusdata);
1041   new(lrecpoints[fNpoint++]) AliPMDrecpoint1(newrecpoint);
1042   delete newrecpoint;
1043 }
1044 // ------------------------------------------------------------------------- //
1045 void AliPMDClusterFinder::AddRecHit(Int_t celldataX,Int_t celldataY,
1046                                     Int_t celldataTr, Int_t celldataPid)
1047 {
1048   // Add associated cell hits to the Reconstructed points
1049   //
1050   TClonesArray &lrechits = *fRechits;
1051   AliPMDrechit *newrechit;
1052   newrechit = new AliPMDrechit(celldataX, celldataY, celldataTr, celldataPid);
1053   new(lrechits[fNhit++]) AliPMDrechit(newrechit);
1054   delete newrechit;
1055 }
1056 // ------------------------------------------------------------------------- //
1057 void AliPMDClusterFinder::ResetCellADC()
1058 {
1059   // Reset the individual cell ADC value to zero
1060   //
1061   for(Int_t irow = 0; irow < fgkRow; irow++)
1062     {
1063       for(Int_t icol = 0; icol < fgkCol; icol++)
1064         {
1065           fCellTrack[irow][icol] = -1;
1066           fCellPid[irow][icol]   = -1;
1067           fCellADC[irow][icol]   = 0.;
1068         }
1069     }
1070 }
1071 // ------------------------------------------------------------------------- //
1072 void AliPMDClusterFinder::ResetRecpoint()
1073 {
1074   // Clear the list of reconstructed points
1075   fNpoint = 0;
1076   if (fRecpoints) fRecpoints->Clear();
1077 }
1078 // ------------------------------------------------------------------------- //
1079 void AliPMDClusterFinder::ResetRechit()
1080 {
1081   // Clear the list of reconstructed points
1082   fNhit = 0;
1083   if (fRechits) fRechits->Clear();
1084 }
1085 // ------------------------------------------------------------------------- //
1086 void AliPMDClusterFinder::Load()
1087 {
1088   // Load all the *.root files
1089   //
1090   fPMDLoader->LoadDigits("READ");
1091   fPMDLoader->LoadRecPoints("recreate");
1092 }
1093 // ------------------------------------------------------------------------- //
1094 void AliPMDClusterFinder::LoadClusters()
1095 {
1096   // Load all the *.root files
1097   //
1098   fPMDLoader->LoadRecPoints("recreate");
1099 }
1100 // ------------------------------------------------------------------------- //
1101 void AliPMDClusterFinder::UnLoad()
1102 {
1103   // Unload all the *.root files
1104   //
1105   fPMDLoader->UnloadDigits();
1106   fPMDLoader->UnloadRecPoints();
1107 }
1108 // ------------------------------------------------------------------------- //
1109 void AliPMDClusterFinder::UnLoadClusters()
1110 {
1111   // Unload all the *.root files
1112   //
1113   fPMDLoader->UnloadRecPoints();
1114 }
1115 // ------------------------------------------------------------------------- //
1116 AliPMDCalibData* AliPMDClusterFinder::GetCalibGain() const
1117 {
1118   // The run number will be centralized in AliCDBManager,
1119   // you don't need to set it here!
1120   // Added by ZA
1121   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Gain");
1122
1123   if(!entry)  AliFatal("Calibration object retrieval failed! ");
1124
1125   AliPMDCalibData *calibdata=0;
1126   if (entry) calibdata = (AliPMDCalibData*) entry->GetObject();
1127
1128   if (!calibdata)  AliFatal("No calibration data from calibration database !");
1129
1130   return calibdata;
1131 }
1132 // ------------------------------------------------------------------------- //
1133 AliPMDPedestal* AliPMDClusterFinder::GetCalibPed() const
1134 {
1135   // The run number will be centralized in AliCDBManager,
1136   // you don't need to set it here!
1137   AliCDBEntry  *entry = AliCDBManager::Instance()->Get("PMD/Calib/Ped");
1138
1139   if(!entry) AliFatal("Pedestal object retrieval failed!");
1140
1141   AliPMDPedestal *pedestal = 0;
1142   if (entry) pedestal = (AliPMDPedestal*) entry->GetObject();
1143
1144   if (!pedestal)  AliFatal("No pedestal data from pedestal database !");
1145
1146   return pedestal;
1147 }