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