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