PMD/AliPMDtracker.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   : March 25 2004                    //
  19 //  This reads the file PMD.RecPoints.root(TreeR),     //
  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 <TMath.h>
  27 #include <TTree.h>
  28 #include <TObjArray.h>
  29 #include <TClonesArray.h>
  30 #include <TFile.h>
  31 #include <TBranch.h>
  32 #include <TNtuple.h>
  33 #include <TParticle.h>
  34
  35 #include "AliPMDcluster.h"
  36 #include "AliPMDclupid.h"
  37 #include "AliPMDrecpoint1.h"
  38 #include "AliPMDrecdata.h"
  39 #include "AliPMDrechit.h"
  40 #include "AliPMDUtility.h"
  41 #include "AliPMDDiscriminator.h"
  42 #include "AliPMDEmpDiscriminator.h"
  43 #include "AliPMDtracker.h"
  44
  45 #include "AliESDPmdTrack.h"
  46 #include "AliESDEvent.h"
  47 #include "AliLog.h"
  48
  49 ClassImp(AliPMDtracker)
  50
  51 AliPMDtracker::AliPMDtracker():
  52   fTreeR(0),
  53   fRecpoints(new TClonesArray("AliPMDrecpoint1", 10)),
  54   fRechits(new TClonesArray("AliPMDrechit", 10)),
  55   fPMDcontin(new TObjArray()),
  56   fPMDcontout(new TObjArray()),
  57   fPMDutil(new AliPMDUtility()),
  58   fPMDrecpoint(0),
  59   fPMDclin(0),
  60   fPMDclout(0),
  61   fXvertex(0.),
  62   fYvertex(0.),
  63   fZvertex(0.),
  64   fSigmaX(0.),
  65   fSigmaY(0.),
  66   fSigmaZ(0.)
  67 {
  68   //
  69   // Default Constructor
  70   //
  71 }
  72 //--------------------------------------------------------------------//
  73 AliPMDtracker:: AliPMDtracker(const AliPMDtracker & /* tracker */):
  74   TObject(/* tracker */),
  75   fTreeR(0),
  76   fRecpoints(NULL),
  77   fRechits(NULL),
  78   fPMDcontin(NULL),
  79   fPMDcontout(NULL),
  80   fPMDutil(NULL),
  81   fPMDrecpoint(0),
  82   fPMDclin(0),
  83   fPMDclout(0),
  84   fXvertex(0.),
  85   fYvertex(0.),
  86   fZvertex(0.),
  87   fSigmaX(0.),
  88   fSigmaY(0.),
  89   fSigmaZ(0.)
  90 {
  91   // copy constructor
  92   AliError("Copy constructor not allowed");
  93 }
  94
  95 //--------------------------------------------------------------------//
  96 AliPMDtracker& AliPMDtracker::operator=(const AliPMDtracker & /* tracker */)
  97 {
  98  // assignment operator
  99   AliError("Assignment operator not allowed");
 100   return *this;
 101 }
 102
 103 //--------------------------------------------------------------------//
 104 AliPMDtracker::~AliPMDtracker()
 105 {
 106   // Destructor
 107   if (fRecpoints)
 108     {
 109       fRecpoints->Clear();
 110     }
 111   if (fRechits)
 112     {
 113       fRechits->Clear();
 114     }
 115
 116   if (fPMDcontin)
 117     {
 118       fPMDcontin->Delete();
 119       delete fPMDcontin;
 120       fPMDcontin=0;
 121
 122     }
 123   if (fPMDcontout)
 124   {
 125       fPMDcontout->Delete();
 126       delete fPMDcontout;
 127       fPMDcontout=0;
 128
 129     }
 130   delete fPMDutil;
 131 }
 132 //--------------------------------------------------------------------//
 133 void AliPMDtracker::LoadClusters(TTree *treein)
 134 {
 135   // Load the Reconstructed tree
 136   fTreeR = treein;
 137 }
 138 //--------------------------------------------------------------------//
 139 void AliPMDtracker::Clusters2Tracks(AliESDEvent *event)
 140 {
 141   // Converts digits to recpoints after running clustering
 142   // algorithm on CPV plane and PREshower plane
 143   //
 144
 145   Int_t   idet;
 146   Int_t   ismn;
 147   Int_t   trackno, trackpid;
 148   Float_t clusdata[6];
 149
 150   Int_t *irow;
 151   Int_t *icol;
 152   Int_t *itra;
 153   Int_t *ipid;
 154   Float_t *cadc;
 155
 156   AliPMDrechit *rechit = 0x0;
 157
 158   TBranch *branch = fTreeR->GetBranch("PMDRecpoint");
 159   if (!branch)
 160     {
 161       AliError("PMDRecpoint branch not found");
 162       return;
 163     }
 164   branch->SetAddress(&fRecpoints);
 165
 166   TBranch *branch1 = fTreeR->GetBranch("PMDRechit");
 167   if (!branch1)
 168     {
 169       AliError("PMDRechit branch not found");
 170       return;
 171     }
 172   branch1->SetAddress(&fRechits);
 173
 174   Int_t   nmodules = (Int_t) branch->GetEntries();
 175
 176   AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
 177   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 178     {
 179       branch->GetEntry(imodule);
 180       Int_t nentries = fRecpoints->GetLast();
 181       AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
 182                       ,nentries));
 183
 184       Int_t ncrhit = 0;
 185
 186       for(Int_t ient = 0; ient < nentries+1; ient++)
 187         {
 188           fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
 189           idet        = fPMDrecpoint->GetDetector();
 190           ismn        = fPMDrecpoint->GetSMNumber();
 191           clusdata[0] = fPMDrecpoint->GetClusX();
 192           clusdata[1] = fPMDrecpoint->GetClusY();
 193           clusdata[2] = fPMDrecpoint->GetClusADC();
 194           clusdata[3] = fPMDrecpoint->GetClusCells();
 195           clusdata[4] = fPMDrecpoint->GetClusSigmaX();
 196           clusdata[5] = fPMDrecpoint->GetClusSigmaY();
 197
 198           if (clusdata[4] != -99. && clusdata[5] != -99.)
 199             {
 200               // extract the associated cell information
 201               branch1->GetEntry(ncrhit);
 202               Int_t nenbr1 = fRechits->GetLast() + 1;
 203
 204               irow = new Int_t[nenbr1];
 205               icol = new Int_t[nenbr1];
 206               itra = new Int_t[nenbr1];
 207               ipid = new Int_t[nenbr1];
 208               cadc = new Float_t[nenbr1];
 209
 210               for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
 211                 {
 212                   rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
 213                   //irow[ient1] = rechit->GetCellX();
 214                   //icol[ient1] = rechit->GetCellY();
 215                   itra[ient1] = rechit->GetCellTrack();
 216                   ipid[ient1] = rechit->GetCellPid();
 217                   cadc[ient1] = rechit->GetCellAdc();
 218                 }
 219               AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
 220                                    trackno, trackpid);
 221
 222               delete [] irow;
 223               delete [] icol;
 224               delete [] itra;
 225               delete [] ipid;
 226               delete [] cadc;
 227
 228               fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
 229               fPMDcontin->Add(fPMDclin);
 230
 231               ncrhit++;
 232             }
 233         }
 234     }
 235
 236   AliPMDDiscriminator *pmddiscriminator = new AliPMDEmpDiscriminator();
 237   pmddiscriminator->Discrimination(fPMDcontin,fPMDcontout);
 238
 239   const Float_t kzpos = 361.5;    // middle of the PMD
 240
 241   Int_t   det,smn,trno,trpid,mstat;
 242   Float_t xpos,ypos;
 243   Float_t adc, ncell, radx, rady;
 244   Float_t xglobal = 0., yglobal = 0., zglobal = 0;
 245   Float_t pid;
 246
 247
 248   Int_t nentries2 = fPMDcontout->GetEntries();
 249   AliDebug(1,Form("Number of clusters coming after discrimination = %d"
 250                   ,nentries2));
 251   for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
 252     {
 253       fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
 254
 255       det   = fPMDclout->GetDetector();
 256       smn   = fPMDclout->GetSMN();
 257       trno  = fPMDclout->GetClusTrackNo();
 258       trpid = fPMDclout->GetClusTrackPid();
 259       mstat = fPMDclout->GetClusMatching();
 260       xpos  = fPMDclout->GetClusX();
 261       ypos  = fPMDclout->GetClusY();
 262       adc   = fPMDclout->GetClusADC();
 263       ncell = fPMDclout->GetClusCells();
 264       radx  = fPMDclout->GetClusSigmaX();
 265       rady  = fPMDclout->GetClusSigmaY();
 266       pid   = fPMDclout->GetClusPID();
 267
 268       //
 269       /**********************************************************************
 270        *    det   : Detector, 0: PRE & 1:CPV                                *
 271        *    smn   : Serial Module Number 0 to 23 for each plane             *
 272        *    xpos  : x-position of the cluster                               *
 273        *    ypos  : y-position of the cluster                               *
 274        *            THESE xpos & ypos are not the true xpos and ypos        *
 275        *            for some of the unit modules. They are rotated.         *
 276        *    adc   : ADC contained in the cluster                            *
 277        *    ncell : Number of cells contained in the cluster                *
 278        *    rad   : radius of the cluster (1d fit)                          *
 279        **********************************************************************/
 280       //
 281
 282       fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal);
 283
 284       if (det == 0)
 285         {
 286           zglobal = kzpos + 1.6; // PREshower plane
 287         }
 288       else if (det == 1)
 289         {
 290           zglobal = kzpos - 1.7; // CPV plane
 291         }
 292
 293       // Fill ESD
 294
 295       AliESDPmdTrack *esdpmdtr = new  AliESDPmdTrack();
 296
 297       esdpmdtr->SetDetector(det);
 298       esdpmdtr->SetSmn(smn);
 299       esdpmdtr->SetClusterTrackNo(trno);
 300       esdpmdtr->SetClusterTrackPid(trpid);
 301       esdpmdtr->SetClusterMatching(mstat);
 302
 303       esdpmdtr->SetClusterX(xglobal);
 304       esdpmdtr->SetClusterY(yglobal);
 305       esdpmdtr->SetClusterZ(zglobal);
 306       esdpmdtr->SetClusterADC(adc);
 307       esdpmdtr->SetClusterCells(ncell);
 308       esdpmdtr->SetClusterPID(pid);
 309       esdpmdtr->SetClusterSigmaX(radx);
 310       esdpmdtr->SetClusterSigmaY(rady);
 311
 312       event->AddPmdTrack(esdpmdtr);
 313     }
 314
 315   fPMDcontin->Delete();
 316   fPMDcontout->Delete();
 317
 318 }
 319 //--------------------------------------------------------------------//
 320 void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
 321                                          Int_t *ipid, Float_t *cadc,
 322                                          Int_t &trackno, Int_t &trackpid)
 323 {
 324   // assign the track number and the corresponding pid to a cluster
 325   // split cluster part will be done at the time of calculating eff/pur
 326
 327   Int_t *phentry = new Int_t [nentry];
 328   Int_t *trenergy = 0x0;
 329   Int_t *sortcoord = 0x0;
 330
 331   Int_t ngtrack = 0;
 332   for (Int_t i = 0; i < nentry; i++)
 333     {
 334       phentry[i] = -1;
 335       if (ipid[i] == 22)
 336         {
 337           phentry[ngtrack] = i;
 338           ngtrack++;
 339         }
 340     }
 341
 342   if (ngtrack == 0)
 343     {
 344       // hadron track
 345       // no need of track number, set to -1
 346       trackpid = 8;
 347       trackno  = -1;
 348     }
 349   else if (ngtrack == 1)
 350     {
 351       // only one photon track
 352       // track number set to photon track
 353       trackpid = 1;
 354       trackno  = itra[phentry[0]];
 355     }
 356   else if (ngtrack > 1)
 357     {
 358       // more than one photon track
 359
 360       trenergy  = new Int_t [ngtrack];
 361       sortcoord = new Int_t [ngtrack];
 362       for (Int_t i = 0; i < ngtrack; i++)
 363         {
 364           trenergy[i] = 0.;
 365           for (Int_t j = 0; j < nentry; j++)
 366             {
 367               if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
 368                 {
 369                   trenergy[i] += cadc[j];
 370                 }
 371             }
 372         }
 373
 374       Bool_t jsort = true;
 375       TMath::Sort(ngtrack,trenergy,sortcoord,jsort);
 376
 377       Int_t gtr = sortcoord[0];
 378       trackno  = itra[phentry[gtr]];   // highest adc track
 379       trackpid = 1;
 380
 381       delete [] trenergy;
 382       delete [] sortcoord;
 383
 384     }   // end of ngtrack > 1
 385
 386
 387
 388 }
 389 //--------------------------------------------------------------------//
 390 void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
 391 {
 392   fXvertex = vtx[0];
 393   fYvertex = vtx[1];
 394   fZvertex = vtx[2];
 395   fSigmaX  = evtx[0];
 396   fSigmaY  = evtx[1];
 397   fSigmaZ  = evtx[2];
 398 }
 399 //--------------------------------------------------------------------//
 400 void AliPMDtracker::ResetClusters()
 401 {
 402   if (fRecpoints) fRecpoints->Clear();
 403 }
 404 //--------------------------------------------------------------------//