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 ncrhit = 0;
 175   Int_t   nmodules = (Int_t) branch->GetEntries();
 176
 177   AliDebug(1,Form("Number of modules filled in treeR = %d",nmodules));
 178   for (Int_t imodule = 0; imodule < nmodules; imodule++)
 179     {
 180       branch->GetEntry(imodule);
 181       Int_t nentries = fRecpoints->GetLast();
 182       AliDebug(2,Form("Number of clusters per modules filled in treeR = %d"
 183                       ,nentries));
 184
 185       for(Int_t ient = 0; ient < nentries+1; ient++)
 186         {
 187           fPMDrecpoint = (AliPMDrecpoint1*)fRecpoints->UncheckedAt(ient);
 188           idet        = fPMDrecpoint->GetDetector();
 189           ismn        = fPMDrecpoint->GetSMNumber();
 190           clusdata[0] = fPMDrecpoint->GetClusX();
 191           clusdata[1] = fPMDrecpoint->GetClusY();
 192           clusdata[2] = fPMDrecpoint->GetClusADC();
 193           clusdata[3] = fPMDrecpoint->GetClusCells();
 194           clusdata[4] = fPMDrecpoint->GetClusSigmaX();
 195           clusdata[5] = fPMDrecpoint->GetClusSigmaY();
 196
 197           if (clusdata[4] < -90. && clusdata[5] < -90.)
 198             {
 199               // extract the associated cell information
 200               branch1->GetEntry(ncrhit);
 201               Int_t nenbr1 = fRechits->GetLast() + 1;
 202
 203               irow = new Int_t[nenbr1];
 204               icol = new Int_t[nenbr1];
 205               itra = new Int_t[nenbr1];
 206               ipid = new Int_t[nenbr1];
 207               cadc = new Float_t[nenbr1];
 208
 209               for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
 210                 {
 211                   rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
 212                   //irow[ient1] = rechit->GetCellX();
 213                   //icol[ient1] = rechit->GetCellY();
 214                   itra[ient1] = rechit->GetCellTrack();
 215                   ipid[ient1] = rechit->GetCellPid();
 216                   cadc[ient1] = rechit->GetCellAdc();
 217                 }
 218               AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
 219                                    trackno, trackpid);
 220
 221               delete [] irow;
 222               delete [] icol;
 223               delete [] itra;
 224               delete [] ipid;
 225               delete [] cadc;
 226
 227               fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
 228               fPMDcontin->Add(fPMDclin);
 229
 230               ncrhit++;
 231             }
 232         }
 233     }
 234
 235   AliPMDDiscriminator *pmddiscriminator = new AliPMDEmpDiscriminator();
 236   pmddiscriminator->Discrimination(fPMDcontin,fPMDcontout);
 237
 238   const Float_t kzpos = 361.5;    // middle of the PMD
 239
 240   Int_t   det,smn,trno,trpid,mstat;
 241   Float_t xpos,ypos;
 242   Float_t adc, ncell, radx, rady;
 243   Float_t xglobal = 0., yglobal = 0., zglobal = 0;
 244   Float_t pid;
 245
 246
 247   Int_t nentries2 = fPMDcontout->GetEntries();
 248   AliDebug(1,Form("Number of clusters coming after discrimination = %d"
 249                   ,nentries2));
 250   for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
 251     {
 252       fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
 253
 254       det   = fPMDclout->GetDetector();
 255       smn   = fPMDclout->GetSMN();
 256       trno  = fPMDclout->GetClusTrackNo();
 257       trpid = fPMDclout->GetClusTrackPid();
 258       mstat = fPMDclout->GetClusMatching();
 259       xpos  = fPMDclout->GetClusX();
 260       ypos  = fPMDclout->GetClusY();
 261       adc   = fPMDclout->GetClusADC();
 262       ncell = fPMDclout->GetClusCells();
 263       radx  = fPMDclout->GetClusSigmaX();
 264       rady  = fPMDclout->GetClusSigmaY();
 265       pid   = fPMDclout->GetClusPID();
 266
 267       //
 268       /**********************************************************************
 269        *    det   : Detector, 0: PRE & 1:CPV                                *
 270        *    smn   : Serial Module Number 0 to 23 for each plane             *
 271        *    xpos  : x-position of the cluster                               *
 272        *    ypos  : y-position of the cluster                               *
 273        *            THESE xpos & ypos are not the true xpos and ypos        *
 274        *            for some of the unit modules. They are rotated.         *
 275        *    adc   : ADC contained in the cluster                            *
 276        *    ncell : Number of cells contained in the cluster                *
 277        *    rad   : radius of the cluster (1d fit)                          *
 278        **********************************************************************/
 279       //
 280
 281       fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal);
 282
 283       if (det == 0)
 284         {
 285           zglobal = kzpos + 1.6; // PREshower plane
 286         }
 287       else if (det == 1)
 288         {
 289           zglobal = kzpos - 1.7; // CPV plane
 290         }
 291
 292       // Fill ESD
 293
 294       AliESDPmdTrack *esdpmdtr = new  AliESDPmdTrack();
 295
 296       esdpmdtr->SetDetector(det);
 297       esdpmdtr->SetSmn(smn);
 298       esdpmdtr->SetClusterTrackNo(trno);
 299       esdpmdtr->SetClusterTrackPid(trpid);
 300       esdpmdtr->SetClusterMatching(mstat);
 301
 302       esdpmdtr->SetClusterX(xglobal);
 303       esdpmdtr->SetClusterY(yglobal);
 304       esdpmdtr->SetClusterZ(zglobal);
 305       esdpmdtr->SetClusterADC(adc);
 306       esdpmdtr->SetClusterCells(ncell);
 307       esdpmdtr->SetClusterPID(pid);
 308       esdpmdtr->SetClusterSigmaX(radx);
 309       esdpmdtr->SetClusterSigmaY(rady);
 310
 311       event->AddPmdTrack(esdpmdtr);
 312     }
 313
 314   fPMDcontin->Delete();
 315   fPMDcontout->Delete();
 316
 317 }
 318 //--------------------------------------------------------------------//
 319 void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
 320                                          Int_t *ipid, Float_t *cadc,
 321                                          Int_t &trackno, Int_t &trackpid)
 322 {
 323   // assign the track number and the corresponding pid to a cluster
 324   // split cluster part will be done at the time of calculating eff/pur
 325
 326   Int_t *phentry = new Int_t [nentry];
 327   Int_t *trenergy = 0x0;
 328   Int_t *sortcoord = 0x0;
 329
 330   Int_t ngtrack = 0;
 331   for (Int_t i = 0; i < nentry; i++)
 332     {
 333       phentry[i] = -1;
 334       if (ipid[i] == 22)
 335         {
 336           phentry[ngtrack] = i;
 337           ngtrack++;
 338         }
 339     }
 340
 341   if (ngtrack == 0)
 342     {
 343       // hadron track
 344       // no need of track number, set to -1
 345       trackpid = 8;
 346       trackno  = -1;
 347     }
 348   else if (ngtrack == 1)
 349     {
 350       // only one photon track
 351       // track number set to photon track
 352       trackpid = 1;
 353       trackno  = itra[phentry[0]];
 354     }
 355   else if (ngtrack > 1)
 356     {
 357       // more than one photon track
 358
 359       trenergy  = new Int_t [ngtrack];
 360       sortcoord = new Int_t [ngtrack];
 361       for (Int_t i = 0; i < ngtrack; i++)
 362         {
 363           trenergy[i] = 0.;
 364           for (Int_t j = 0; j < nentry; j++)
 365             {
 366               if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
 367                 {
 368                   trenergy[i] += cadc[j];
 369                 }
 370             }
 371         }
 372
 373       Bool_t jsort = true;
 374       TMath::Sort(ngtrack,trenergy,sortcoord,jsort);
 375
 376       Int_t gtr = sortcoord[0];
 377       trackno  = itra[phentry[gtr]];   // highest adc track
 378       trackpid = 1;
 379
 380       delete [] trenergy;
 381       delete [] sortcoord;
 382
 383     }   // end of ngtrack > 1
 384
 385
 386
 387 }
 388 //--------------------------------------------------------------------//
 389 void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
 390 {
 391   fXvertex = vtx[0];
 392   fYvertex = vtx[1];
 393   fZvertex = vtx[2];
 394   fSigmaX  = evtx[0];
 395   fSigmaY  = evtx[1];
 396   fSigmaZ  = evtx[2];
 397 }
 398 //--------------------------------------------------------------------//
 399 void AliPMDtracker::ResetClusters()
 400 {
 401   if (fRecpoints) fRecpoints->Clear();
 402 }
 403 //--------------------------------------------------------------------//