PMD/AliPMDRecpointRead.C

   1 // ----------------------------------------------------//
   2 //                                                     //
   3 //       This macro reads the PMD clusters which       //
   4 //       are stored in the file "PMD.RecPoints.root"   //
   5 //                                                     //
   6 // ----------------------------------------------------//
   7
   8 #include <Riostream.h>
   9 #include "TBranch.h"
  10 #include "TStopwatch.h"
  11
  12 extern AliRun *gAlice;
  13
  14 Int_t AliPMDRecpointRead(Int_t nevent = 1)
  15 {
  16   if (gAlice)
  17     {
  18       delete gAlice->GetRunLoader();
  19       delete gAlice;//if everything was OK here it is already NULL
  20       gAlice = 0x0;
  21     }
  22   AliRunLoader *fRunLoader = AliRunLoader::Open("galice.root","Event","update");
  23   if (!fRunLoader)
  24     {
  25       cerr<<"Can't load RunLoader"<<endl;
  26       return 1;
  27     }
  28   AliLoader *pmdloader = fRunLoader->GetLoader("PMDLoader");
  29   Int_t nevent = fRunLoader->GetNumberOfEvents();
  30   cout << " * *********** nevent = " << nevent << endl;
  31
  32   if (pmdloader == 0x0)
  33     {
  34       cerr<<" ===> Can not find PMD or PMDLoader <===\n";
  35       delete fRunLoader;
  36       return 2;
  37     }
  38
  39   pmdloader->LoadRecPoints("READ");
  40   TClonesArray *fRecpoints;
  41   AliPMDUtility *cc = new AliPMDUtility();
  42   TH2F *h2 = new TH2F("h2"," ",100,-100.,100.,100,-100.,100.);
  43
  44   FILE *fpw = fopen("junk_rec.dat","w");
  45
  46   for (Int_t ievt = 0; ievt < nevent; ievt++)
  47     {
  48       fRunLoader->GetEvent(ievt);
  49       TTree *treeR = pmdloader->TreeR();
  50       if (treeR == 0x0)
  51         {
  52           cout << " Can not get TreeR" << endl;
  53           return 3;
  54         }
  55
  56       AliPMDrecpoint1  *pmdrecpoint;
  57       TBranch *branch = treeR->GetBranch("PMDRecpoint");
  58       branch->SetAddress(&fRecpoints);
  59       /**********************************************************************
  60        *    det   : Detector, 0: PRE & 1:CPV                                *
  61        *    smn   : Serial Module Number from which Super Module Number     *
  62        *            and Unit Module Numbers are extracted                   *
  63        *    xpos  : x-position of the cluster                               *
  64        *    ypos  : y-position of the cluster                               *
  65        *            THESE xpos & ypos are not the true xpos and ypos        *
  66        *            for some of the unit modules. They are rotated.         *
  67        *    adc   : ADC contained in the cluster                            *
  68        *    ncell : Number of cells contained in the cluster                *
  69        *    rad   : radius of the cluster (1d fit)                          *
  70        *    ism   : Supermodule number extracted from smn                   *
  71        *    ium   : Unit module number extracted from smn                   *
  72        *    xpad  : TRUE x-position of the cluster                          *
  73        *    ypad  : TRUE y-position of the cluster                          *
  74        **********************************************************************/
  75       Int_t   ism, ium;
  76       Int_t   det,smn;
  77       Float_t xpos,ypos, xpad, ypad;
  78       Float_t adc, ncell, rad;
  79       Float_t xx, yy;
  80       Int_t   nmodules = treeR->GetEntries();
  81       cout << " nmodules = " << nmodules << endl;
  82       for (Int_t imodule = 0; imodule < nmodules; imodule++)
  83         {
  84           treeR->GetEntry(imodule);
  85           Int_t nentries = fRecpoints->GetLast();
  86           for(Int_t ient = 0; ient < nentries+1; ient++)
  87             {
  88               pmdrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
  89               det   = (Int_t) pmdrecpoint->GetDetector();
  90               smn   = (Int_t) pmdrecpoint->GetSMNumber();
  91               xpos  = pmdrecpoint->GetClusX();
  92               ypos  = pmdrecpoint->GetClusY();
  93               adc   = pmdrecpoint->GetClusADC();
  94               ncell = pmdrecpoint->GetClusCells();
  95               rad   = pmdrecpoint->GetClusRadius();
  96               //
  97               // Now change the xpos and ypos to its original values
  98               // for the unit modules which are earlier changed.
  99               // xpad and ypad are the real positions.
 100               //
 101               if(det == 0 || det == 1)
 102                 {
 103                   if(smn < 12)
 104                     {
 105                       ism  = smn/6;
 106                       ium  = smn - ism*6;
 107                       xpad = ypos;
 108                       ypad = xpos;
 109                     }
 110                   else if( smn >= 12 && smn < 24)
 111                     {
 112                       ism  = smn/6;
 113                       ium  = smn - ism*6;
 114                       xpad = xpos;
 115                       ypad = ypos;
 116                     }
 117                 }
 118               //
 119               // User has to plug in his analysis code here
 120               //
 121
 122               fprintf(fpw,"%d %d %d %d %f %f %f %f %f\n",
 123                       det,smn,ism,ium,xpad,ypad,adc,ncell,rad);
 124               //
 125               // Plot the cluster centroid to see the PMD geometry
 126               // using the PMD Utility class
 127               //
 128               if (det == 1)
 129                 {
 130                   // Draw only for PRE plane
 131                   cc->RectGeomCellPos(ism,ium,xpad,ypad,xx,yy);
 132                   h2->Fill(xx,yy);
 133                 }
 134
 135               //
 136               // End of the User code
 137               //
 138             }
 139         }
 140
 141     }
 142
 143   h2->Draw();
 144
 145   fclose(fpw);
 146   return 0;
 147 }
 148