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