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 "AliPMDcluster.h"
  39 #include "AliPMDrecpoint1.h"
  40 #include "AliPMDRawStream.h"
  41
  42
  43
  44 ClassImp(AliPMDClusterFinder)
  45
  46 AliPMDClusterFinder::AliPMDClusterFinder():
  47   fRunLoader(0),
  48   fPMDLoader(0),
  49   fTreeD(0),
  50   fTreeR(0),
  51   fDigits(new TClonesArray("AliPMDdigit", 1000)),
  52   fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
  53   fNpoint(0),
  54   fEcut(0.)
  55 {
  56 //
  57 // Constructor
  58 //
  59 }
  60 // ------------------------------------------------------------------------- //
  61 AliPMDClusterFinder::AliPMDClusterFinder(AliRunLoader* runLoader):
  62   fRunLoader(runLoader),
  63   fPMDLoader(runLoader->GetLoader("PMDLoader")),
  64   fTreeD(0),
  65   fTreeR(0),
  66   fDigits(new TClonesArray("AliPMDdigit", 1000)),
  67   fRecpoints(new TClonesArray("AliPMDrecpoint1", 1000)),
  68   fNpoint(0),
  69   fEcut(0.)
  70 {
  71 //
  72 // Constructor
  73 //
  74 }
  75 // ------------------------------------------------------------------------- //
  76 AliPMDClusterFinder::~AliPMDClusterFinder()
  77 {
  78   // Destructor
  79   if (fDigits)
  80     {
  81       fDigits->Delete();
  82       delete fDigits;
  83       fDigits=0;
  84     }
  85   if (fRecpoints)
  86     {
  87       fRecpoints->Delete();
  88       delete fRecpoints;
  89       fRecpoints=0;
  90     }
  91 }
  92 // ------------------------------------------------------------------------- //
  93
  94 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt)
  95 {
  96   // Converts digits to recpoints after running clustering
  97   // algorithm on CPV plane and PREshower plane
  98   //
  99   Int_t    det = 0,smn = 0;
 100   Int_t    xpos,ypos;
 101   Float_t  adc;
 102   Int_t    ismn;
 103   Int_t    idet;
 104   Float_t  clusdata[5];
 105
 106   TObjArray *pmdcont = new TObjArray();
 107   AliPMDClustering *pmdclust = new AliPMDClustering();
 108
 109   pmdclust->SetEdepCut(fEcut);
 110
 111   fRunLoader->GetEvent(ievt);
 112
 113
 114   fTreeD = fPMDLoader->TreeD();
 115   if (fTreeD == 0x0)
 116     {
 117       AliFatal("AliPMDClusterFinder: Can not get TreeD");
 118
 119     }
 120   AliPMDdigit  *pmddigit;
 121   TBranch *branch = fTreeD->GetBranch("PMDDigit");
 122   branch->SetAddress(&fDigits);
 123
 124   ResetRecpoint();
 125
 126   fTreeR = fPMDLoader->TreeR();
 127   if (fTreeR == 0x0)
 128     {
 129       fPMDLoader->MakeTree("R");
 130       fTreeR = fPMDLoader->TreeR();
 131     }
 132
 133   Int_t bufsize = 16000;
 134   fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize);
 135
 136   Int_t nmodules = (Int_t) fTreeD->GetEntries();
 137
 138   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 139     {
 140       ResetCellADC();
 141       fTreeD->GetEntry(imodule);
 142       Int_t nentries = fDigits->GetLast();
 143       for (Int_t ient = 0; ient < nentries+1; ient++)
 144         {
 145           pmddigit = (AliPMDdigit*)fDigits->UncheckedAt(ient);
 146
 147           det    = pmddigit->GetDetector();
 148           smn    = pmddigit->GetSMNumber();
 149           xpos   = pmddigit->GetRow();
 150           ypos   = pmddigit->GetColumn();
 151           adc    = pmddigit->GetADC();
 152           //Int_t trno   = pmddigit->GetTrackNumber();
 153           fCellADC[xpos][ypos] = (Double_t) adc;
 154         }
 155
 156       idet = det;
 157       ismn = smn;
 158       pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
 159
 160       Int_t nentries1 = pmdcont->GetEntries();
 161
 162       AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 163
 164       for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 165         {
 166           AliPMDcluster *pmdcl = (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
 167           idet        = pmdcl->GetDetector();
 168           ismn        = pmdcl->GetSMN();
 169           clusdata[0] = pmdcl->GetClusX();
 170           clusdata[1] = pmdcl->GetClusY();
 171           clusdata[2] = pmdcl->GetClusADC();
 172           clusdata[3] = pmdcl->GetClusCells();
 173           clusdata[4] = pmdcl->GetClusRadius();
 174
 175           AddRecPoint(idet,ismn,clusdata);
 176         }
 177       pmdcont->Clear();
 178
 179       fTreeR->Fill();
 180       ResetRecpoint();
 181
 182     } // modules
 183
 184   ResetCellADC();
 185   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 186   fPMDLoader->WriteRecPoints("OVERWRITE");
 187
 188   //   delete the pointers
 189   delete pmdclust;
 190   delete pmdcont;
 191
 192 }
 193 // ------------------------------------------------------------------------- //
 194
 195 void AliPMDClusterFinder::Digits2RecPoints(AliRawReader *rawReader,
 196                                            TTree *clustersTree)
 197 {
 198   // Converts RAW data to recpoints after running clustering
 199   // algorithm on CPV and PREshower plane
 200   //
 201
 202   Float_t  clusdata[5];
 203
 204   TObjArray *pmdcont = new TObjArray();
 205   AliPMDClustering *pmdclust = new AliPMDClustering();
 206
 207   pmdclust->SetEdepCut(fEcut);
 208
 209   ResetRecpoint();
 210
 211   Int_t bufsize = 16000;
 212   clustersTree->Branch("PMDRecpoint", &fRecpoints, bufsize);
 213
 214   const Int_t kDDL = 6;
 215   const Int_t kRow = 48;
 216   const Int_t kCol = 96;
 217
 218   Int_t idet = 0;
 219   Int_t iSMN = 0;
 220
 221   for (Int_t indexDDL = 0; indexDDL < kDDL; indexDDL++)
 222     {
 223       if (indexDDL < 4)
 224         {
 225           iSMN = 6;
 226         }
 227       else if (indexDDL >= 4)
 228         {
 229           iSMN = 12;
 230         }
 231       Int_t ***precpvADC;
 232       precpvADC = new int **[iSMN];
 233       for (Int_t i=0; i<iSMN; i++) precpvADC[i] = new int *[kRow];
 234       for (Int_t i=0; i<iSMN;i++)
 235         {
 236           for (Int_t j=0; j<kRow; j++) precpvADC[i][j] = new int [kCol];
 237         }
 238       for (Int_t i = 0; i < iSMN; i++)
 239         {
 240           for (Int_t j = 0; j < kRow; j++)
 241             {
 242               for (Int_t k = 0; k < kCol; k++)
 243                 {
 244                   precpvADC[i][j][k] = 0;
 245                 }
 246             }
 247         }
 248       ResetCellADC();
 249       rawReader->Reset();
 250       AliPMDRawStream pmdinput(rawReader);
 251       rawReader->Select(12, indexDDL, indexDDL);
 252       while(pmdinput.Next())
 253         {
 254           Int_t det = pmdinput.GetDetector();
 255           Int_t smn = pmdinput.GetSMN();
 256           //Int_t mcm = pmdinput.GetMCM();
 257           //Int_t chno = pmdinput.GetChannel();
 258           Int_t row = pmdinput.GetRow();
 259           Int_t col = pmdinput.GetColumn();
 260           Int_t sig = pmdinput.GetSignal();
 261
 262           Int_t indexsmn = 0;
 263
 264           if (indexDDL < 4)
 265             {
 266               if (det != 0)
 267                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 268                               indexDDL, det));
 269               indexsmn = smn - indexDDL * 6;
 270             }
 271           else if (indexDDL == 4)
 272             {
 273               if (det != 1)
 274                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 275                               indexDDL, det));
 276               if (smn < 6)
 277                 {
 278                   indexsmn = smn;
 279                 }
 280               else if (smn >= 12 && smn < 18)
 281                 {
 282                   indexsmn = smn - 6;
 283                 }
 284             }
 285           else if (indexDDL == 5)
 286             {
 287               if (det != 1)
 288                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 289                               indexDDL, det));
 290               if (smn >= 6 && smn < 12)
 291                 {
 292                   indexsmn = smn - 6;
 293                 }
 294               else if (smn >= 18 && smn < 24)
 295                 {
 296                   indexsmn = smn - 12;
 297                 }
 298             }
 299           precpvADC[indexsmn][row][col] = sig;
 300         } // while loop
 301
 302       Int_t ismn = 0;
 303       for (Int_t indexsmn = 0; indexsmn < iSMN; indexsmn++)
 304         {
 305           ResetCellADC();
 306           for (Int_t irow = 0; irow < kRow; irow++)
 307             {
 308               for (Int_t icol = 0; icol < kCol; icol++)
 309                 {
 310                   fCellADC[irow][icol] =
 311                     (Double_t) precpvADC[indexsmn][irow][icol];
 312                 } // row
 313             }     // col
 314           if (indexDDL < 4)
 315             {
 316               ismn = indexsmn + indexDDL * 6;
 317               idet = 0;
 318             }
 319           else if (indexDDL == 4)
 320             {
 321               if (indexsmn < 6)
 322                 {
 323                   ismn = indexsmn;
 324                 }
 325               else if (indexsmn >= 6 && indexsmn < 12)
 326                 {
 327                   ismn = indexsmn + 6;
 328                 }
 329               idet = 1;
 330             }
 331           else if (indexDDL == 5)
 332             {
 333               if (indexsmn < 6)
 334                 {
 335                   ismn = indexsmn + 6;
 336                 }
 337               else if (indexsmn >= 6 && indexsmn < 12)
 338                 {
 339                   ismn = indexsmn + 12;
 340                 }
 341               idet = 1;
 342             }
 343
 344
 345           pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
 346           Int_t nentries1 = pmdcont->GetEntries();
 347
 348           AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 349
 350           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 351             {
 352               AliPMDcluster *pmdcl =
 353                 (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
 354               idet        = pmdcl->GetDetector();
 355               ismn        = pmdcl->GetSMN();
 356               clusdata[0] = pmdcl->GetClusX();
 357               clusdata[1] = pmdcl->GetClusY();
 358               clusdata[2] = pmdcl->GetClusADC();
 359               clusdata[3] = pmdcl->GetClusCells();
 360               clusdata[4] = pmdcl->GetClusRadius();
 361
 362               AddRecPoint(idet,ismn,clusdata);
 363             }
 364           pmdcont->Clear();
 365
 366           //fTreeR->Fill();
 367           clustersTree->Fill();
 368           ResetRecpoint();
 369
 370
 371         } // smn
 372
 373       for (Int_t i=0; i<iSMN; i++)
 374         {
 375           for (Int_t j=0; j<kRow; j++) delete [] precpvADC[i][j];
 376         }
 377       for (Int_t i=0; i<iSMN; i++) delete [] precpvADC[i];
 378       delete precpvADC;
 379     } // DDL Loop
 380
 381   ResetCellADC();
 382
 383   //   delete the pointers
 384   delete pmdclust;
 385   delete pmdcont;
 386
 387 }
 388 // ------------------------------------------------------------------------- //
 389
 390 void AliPMDClusterFinder::Digits2RecPoints(Int_t ievt, AliRawReader *rawReader)
 391 {
 392   // Converts RAW data to recpoints after running clustering
 393   // algorithm on CPV and PREshower plane
 394   //
 395
 396   Float_t  clusdata[5];
 397
 398   TObjArray *pmdcont = new TObjArray();
 399   AliPMDClustering *pmdclust = new AliPMDClustering();
 400
 401   pmdclust->SetEdepCut(fEcut);
 402
 403   fRunLoader->GetEvent(ievt);
 404
 405   ResetRecpoint();
 406
 407   fTreeR = fPMDLoader->TreeR();
 408   if (fTreeR == 0x0)
 409     {
 410       fPMDLoader->MakeTree("R");
 411       fTreeR = fPMDLoader->TreeR();
 412     }
 413   Int_t bufsize = 16000;
 414   fTreeR->Branch("PMDRecpoint", &fRecpoints, bufsize);
 415
 416   const Int_t kDDL = 6;
 417   const Int_t kRow = 48;
 418   const Int_t kCol = 96;
 419
 420   Int_t idet = 0;
 421   Int_t iSMN = 0;
 422
 423   for (Int_t indexDDL = 0; indexDDL < kDDL; indexDDL++)
 424     {
 425       if (indexDDL < 4)
 426         {
 427           iSMN = 6;
 428         }
 429       else if (indexDDL >= 4)
 430         {
 431           iSMN = 12;
 432         }
 433       Int_t ***precpvADC;
 434       precpvADC = new int **[iSMN];
 435       for (Int_t i=0; i<iSMN; i++) precpvADC[i] = new int *[kRow];
 436       for (Int_t i=0; i<iSMN;i++)
 437         {
 438           for (Int_t j=0; j<kRow; j++) precpvADC[i][j] = new int [kCol];
 439         }
 440       for (Int_t i = 0; i < iSMN; i++)
 441         {
 442           for (Int_t j = 0; j < kRow; j++)
 443             {
 444               for (Int_t k = 0; k < kCol; k++)
 445                 {
 446                   precpvADC[i][j][k] = 0;
 447                 }
 448             }
 449         }
 450       ResetCellADC();
 451       rawReader->Reset();
 452       AliPMDRawStream pmdinput(rawReader);
 453       rawReader->Select(12, indexDDL, indexDDL);
 454       while(pmdinput.Next())
 455         {
 456           Int_t det = pmdinput.GetDetector();
 457           Int_t smn = pmdinput.GetSMN();
 458           //Int_t mcm = pmdinput.GetMCM();
 459           //Int_t chno = pmdinput.GetChannel();
 460           Int_t row = pmdinput.GetRow();
 461           Int_t col = pmdinput.GetColumn();
 462           Int_t sig = pmdinput.GetSignal();
 463
 464           Int_t indexsmn = 0;
 465
 466           if (indexDDL < 4)
 467             {
 468               if (det != 0)
 469                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 470                               indexDDL, det));
 471               indexsmn = smn - indexDDL * 6;
 472             }
 473           else if (indexDDL == 4)
 474             {
 475               if (det != 1)
 476                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 477                               indexDDL, det));
 478               if (smn < 6)
 479                 {
 480                   indexsmn = smn;
 481                 }
 482               else if (smn >= 12 && smn < 18)
 483                 {
 484                   indexsmn = smn - 6;
 485                 }
 486             }
 487           else if (indexDDL == 5)
 488             {
 489               if (det != 1)
 490                 AliError(Form("*DDL %d and Detector NUMBER %d NOT MATCHING *",
 491                               indexDDL, det));
 492               if (smn >= 6 && smn < 12)
 493                 {
 494                   indexsmn = smn - 6;
 495                 }
 496               else if (smn >= 18 && smn < 24)
 497                 {
 498                   indexsmn = smn - 12;
 499                 }
 500             }
 501           precpvADC[indexsmn][row][col] = sig;
 502         } // while loop
 503
 504       Int_t ismn = 0;
 505       for (Int_t indexsmn = 0; indexsmn < iSMN; indexsmn++)
 506         {
 507           ResetCellADC();
 508           for (Int_t irow = 0; irow < kRow; irow++)
 509             {
 510               for (Int_t icol = 0; icol < kCol; icol++)
 511                 {
 512                   fCellADC[irow][icol] =
 513                     (Double_t) precpvADC[indexsmn][irow][icol];
 514                 } // row
 515             }     // col
 516           if (indexDDL < 4)
 517             {
 518               ismn = indexsmn + indexDDL * 6;
 519               idet = 0;
 520             }
 521           else if (indexDDL == 4)
 522             {
 523               if (indexsmn < 6)
 524                 {
 525                   ismn = indexsmn;
 526                 }
 527               else if (indexsmn >= 6 && indexsmn < 12)
 528                 {
 529                   ismn = indexsmn + 6;
 530                 }
 531               idet = 1;
 532             }
 533           else if (indexDDL == 5)
 534             {
 535               if (indexsmn < 6)
 536                 {
 537                   ismn = indexsmn + 6;
 538                 }
 539               else if (indexsmn >= 6 && indexsmn < 12)
 540                 {
 541                   ismn = indexsmn + 12;
 542                 }
 543               idet = 1;
 544             }
 545
 546
 547           pmdclust->DoClust(idet,ismn,fCellADC,pmdcont);
 548           Int_t nentries1 = pmdcont->GetEntries();
 549
 550           AliDebug(1,Form("Total number of clusters/module = %d",nentries1));
 551
 552           for (Int_t ient1 = 0; ient1 < nentries1; ient1++)
 553             {
 554               AliPMDcluster *pmdcl =
 555                 (AliPMDcluster*)pmdcont->UncheckedAt(ient1);
 556               idet        = pmdcl->GetDetector();
 557               ismn        = pmdcl->GetSMN();
 558               clusdata[0] = pmdcl->GetClusX();
 559               clusdata[1] = pmdcl->GetClusY();
 560               clusdata[2] = pmdcl->GetClusADC();
 561               clusdata[3] = pmdcl->GetClusCells();
 562               clusdata[4] = pmdcl->GetClusRadius();
 563
 564               AddRecPoint(idet,ismn,clusdata);
 565             }
 566           pmdcont->Clear();
 567
 568           fTreeR->Fill();
 569           ResetRecpoint();
 570
 571
 572         } // smn
 573
 574       for (Int_t i=0; i<iSMN; i++)
 575         {
 576           for (Int_t j=0; j<kRow; j++) delete [] precpvADC[i][j];
 577         }
 578       for (Int_t i=0; i<iSMN; i++) delete [] precpvADC[i];
 579       delete precpvADC;
 580     } // DDL Loop
 581
 582   ResetCellADC();
 583
 584   fPMDLoader = fRunLoader->GetLoader("PMDLoader");
 585   fPMDLoader->WriteRecPoints("OVERWRITE");
 586
 587   //   delete the pointers
 588   delete pmdclust;
 589   delete pmdcont;
 590
 591 }
 592 // ------------------------------------------------------------------------- //
 593 void AliPMDClusterFinder::SetCellEdepCut(Float_t ecut)
 594 {
 595   fEcut = ecut;
 596 }
 597 // ------------------------------------------------------------------------- //
 598 void AliPMDClusterFinder::AddRecPoint(Int_t idet,Int_t ismn,Float_t *clusdata)
 599 {
 600   // Add Reconstructed points
 601   //
 602   TClonesArray &lrecpoints = *fRecpoints;
 603   AliPMDrecpoint1 *newrecpoint;
 604   newrecpoint = new AliPMDrecpoint1(idet, ismn, clusdata);
 605   new(lrecpoints[fNpoint++]) AliPMDrecpoint1(newrecpoint);
 606   delete newrecpoint;
 607 }
 608 // ------------------------------------------------------------------------- //
 609 void AliPMDClusterFinder::ResetCellADC()
 610 {
 611   // Reset the individual cell ADC value to zero
 612   //
 613   for(Int_t irow = 0; irow < fgkRow; irow++)
 614     {
 615       for(Int_t icol = 0; icol < fgkCol; icol++)
 616         {
 617           fCellADC[irow][icol] = 0.;
 618         }
 619     }
 620 }
 621 // ------------------------------------------------------------------------- //
 622
 623 void AliPMDClusterFinder::ResetRecpoint()
 624 {
 625   // Clear the list of reconstructed points
 626   fNpoint = 0;
 627   if (fRecpoints) fRecpoints->Clear();
 628 }
 629 // ------------------------------------------------------------------------- //
 630 void AliPMDClusterFinder::Load()
 631 {
 632   // Load all the *.root files
 633   //
 634   fPMDLoader->LoadDigits("READ");
 635   fPMDLoader->LoadRecPoints("recreate");
 636 }
 637 // ------------------------------------------------------------------------- //
 638 void AliPMDClusterFinder::LoadClusters()
 639 {
 640   // Load all the *.root files
 641   //
 642   fPMDLoader->LoadRecPoints("recreate");
 643 }
 644 // ------------------------------------------------------------------------- //
 645 void AliPMDClusterFinder::UnLoad()
 646 {
 647   // Unload all the *.root files
 648   //
 649   fPMDLoader->UnloadDigits();
 650   fPMDLoader->UnloadRecPoints();
 651 }
 652 // ------------------------------------------------------------------------- //
 653 void AliPMDClusterFinder::UnLoadClusters()
 654 {
 655   // Unload all the *.root files
 656   //
 657   fPMDLoader->UnloadRecPoints();
 658 }
 659 // ------------------------------------------------------------------------- //