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                   irow[ient1] = -99;
 216                   icol[ient1] = -99;
 217                   itra[ient1] = -99;
 218                   ipid[ient1] = -99;
 219                   cadc[ient1] = 0.;
 220                 }
 221               for (Int_t ient1 = 0; ient1 < nenbr1; ient1++)
 222                 {
 223                   rechit = (AliPMDrechit*)fRechits->UncheckedAt(ient1);
 224                   //irow[ient1] = rechit->GetCellX();
 225                   //icol[ient1] = rechit->GetCellY();
 226                   itra[ient1] = rechit->GetCellTrack();
 227                   ipid[ient1] = rechit->GetCellPid();
 228                   cadc[ient1] = rechit->GetCellAdc();
 229                 }
 230               if (idet == 0)
 231                 {
 232                   AssignTrPidToCluster(nenbr1, itra, ipid, cadc,
 233                                        trackno, trackpid);
 234                 }
 235               else if (idet == 1)
 236                 {
 237                   trackno  = itra[0];
 238                   trackpid = ipid[0];
 239                 }
 240
 241               delete [] irow;
 242               delete [] icol;
 243               delete [] itra;
 244               delete [] ipid;
 245               delete [] cadc;
 246
 247               fPMDclin = new AliPMDrecdata(idet,ismn,trackno,trackpid,clusdata);
 248               fPMDcontin->Add(fPMDclin);
 249
 250               ncrhit++;
 251             }
 252         }
 253     }
 254
 255   AliPMDEmpDiscriminator pmddiscriminator;
 256   pmddiscriminator.Discrimination(fPMDcontin,fPMDcontout);
 257
 258   // alignment implemention
 259
 260   Double_t sectr[4][3] = { {0.,0.,0.},{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
 261   TString snsector="PMD/Sector";
 262   TString symname;
 263   TGeoHMatrix gpmdor;
 264
 265   for(Int_t isector=1; isector<=4; isector++)
 266     {
 267       symname = snsector;
 268       symname += isector;
 269       TGeoHMatrix *gpmdal = AliGeomManager::GetMatrix(symname);
 270       Double_t *tral = gpmdal->GetTranslation();
 271
 272       AliGeomManager::GetOrigGlobalMatrix(symname, gpmdor);
 273       Double_t *tror = gpmdor.GetTranslation();
 274
 275       for(Int_t ixyz=0; ixyz<3; ixyz++)
 276         {
 277           sectr[isector-1][ixyz] = tral[ixyz] - tror[ixyz];
 278         }
 279     }
 280
 281   const Float_t kzpos = 361.5;    // middle of the PMD
 282
 283   Int_t   ix = -1, iy = -1;
 284   Int_t   det = 0, smn = 0, trno = 1, trpid = 0, mstat = 0;
 285   Float_t xpos = 0., ypos = 0.;
 286   Float_t adc = 0., ncell = 0., radx = 0., rady = 0.;
 287   Float_t xglobal = 0., yglobal = 0., zglobal = 0;
 288   Float_t pid = 0.;
 289
 290   fPMDutil->ApplyAlignment(sectr);
 291
 292   Int_t nentries2 = fPMDcontout->GetEntries();
 293   AliDebug(1,Form("Number of clusters coming after discrimination = %d"
 294                   ,nentries2));
 295   for (Int_t ient1 = 0; ient1 < nentries2; ient1++)
 296     {
 297       fPMDclout = (AliPMDclupid*)fPMDcontout->UncheckedAt(ient1);
 298
 299       det   = fPMDclout->GetDetector();
 300       smn   = fPMDclout->GetSMN();
 301       trno  = fPMDclout->GetClusTrackNo();
 302       trpid = fPMDclout->GetClusTrackPid();
 303       mstat = fPMDclout->GetClusMatching();
 304       xpos  = fPMDclout->GetClusX();
 305       ypos  = fPMDclout->GetClusY();
 306       adc   = fPMDclout->GetClusADC();
 307       ncell = fPMDclout->GetClusCells();
 308       radx  = fPMDclout->GetClusSigmaX();
 309       // Here in the variable "rady" we are keeping the row and col
 310       // of the single isolated cluster having ncell = 1 for offline
 311       // calibration
 312
 313       if ((radx > 999. && radx < 1000.) && ncell == 1)
 314         {
 315           if (smn < 12)
 316             {
 317               ix = (Int_t) (ypos +0.5);
 318               iy = (Int_t) xpos;
 319             }
 320           else if (smn > 12 && smn < 24)
 321             {
 322               ix = (Int_t) xpos;
 323               iy = (Int_t) (ypos +0.5);
 324             }
 325           rady = (Float_t) (ix*100 + iy);
 326         }
 327       else
 328         {
 329           rady  = fPMDclout->GetClusSigmaY();
 330         }
 331       pid   = fPMDclout->GetClusPID();
 332
 333       //
 334       /**********************************************************************
 335        *    det   : Detector, 0: PRE & 1:CPV                                *
 336        *    smn   : Serial Module Number 0 to 23 for each plane             *
 337        *    xpos  : x-position of the cluster                               *
 338        *    ypos  : y-position of the cluster                               *
 339        *            THESE xpos & ypos are not the true xpos and ypos        *
 340        *            for some of the unit modules. They are rotated.         *
 341        *    adc   : ADC contained in the cluster                            *
 342        *    ncell : Number of cells contained in the cluster                *
 343        *    rad   : radius of the cluster (1d fit)                          *
 344        **********************************************************************/
 345       //
 346
 347
 348       if (det == 0)
 349         {
 350           zglobal = kzpos + 1.65; // PREshower plane
 351         }
 352       else if (det == 1)
 353         {
 354           zglobal = kzpos - 1.65; // CPV plane
 355         }
 356
 357       fPMDutil->RectGeomCellPos(smn,xpos,ypos,xglobal,yglobal,zglobal);
 358
 359       // Fill ESD
 360
 361       AliESDPmdTrack *esdpmdtr = new  AliESDPmdTrack();
 362
 363       esdpmdtr->SetDetector(det);
 364       esdpmdtr->SetSmn(smn);
 365       esdpmdtr->SetClusterTrackNo(trno);
 366       esdpmdtr->SetClusterTrackPid(trpid);
 367       esdpmdtr->SetClusterMatching(mstat);
 368
 369       esdpmdtr->SetClusterX(xglobal);
 370       esdpmdtr->SetClusterY(yglobal);
 371       esdpmdtr->SetClusterZ(zglobal);
 372       esdpmdtr->SetClusterADC(adc);
 373       esdpmdtr->SetClusterCells(ncell);
 374       esdpmdtr->SetClusterPID(pid);
 375       esdpmdtr->SetClusterSigmaX(radx);
 376       esdpmdtr->SetClusterSigmaY(rady);
 377
 378       event->AddPmdTrack(esdpmdtr);
 379       delete esdpmdtr;
 380     }
 381
 382   fPMDcontin->Delete();
 383   fPMDcontout->Delete();
 384
 385 }
 386 //--------------------------------------------------------------------//
 387 void AliPMDtracker::AssignTrPidToCluster(Int_t nentry, Int_t *itra,
 388                                          Int_t *ipid, Float_t *cadc,
 389                                          Int_t &trackno, Int_t &trackpid)
 390 {
 391   // assign the track number and the corresponding pid to a cluster
 392   // split cluster part will be done at the time of calculating eff/pur
 393
 394   Int_t *phentry = new Int_t [nentry];
 395   Int_t *hadentry = new Int_t [nentry];
 396   Int_t *trpid    = 0x0;
 397   Int_t *sortcoord = 0x0;
 398   Float_t *trenergy = 0x0;
 399
 400   Int_t ngtrack = 0;
 401   Int_t nhtrack = 0;
 402   for (Int_t i = 0; i < nentry; i++)
 403     {
 404       phentry[i] = -1;
 405       hadentry[i] = -1;
 406
 407       if (ipid[i] == 22)
 408         {
 409           phentry[ngtrack] = i;
 410           ngtrack++;
 411         }
 412       else if (ipid[i] != 22)
 413         {
 414           hadentry[nhtrack] = i;
 415           nhtrack++;
 416         }
 417     }
 418
 419   Int_t nghadtrack = ngtrack + nhtrack;
 420
 421   if (ngtrack == 0)
 422     {
 423       // hadron track
 424       // no need of track number, set to -1
 425       trackpid = 8;
 426       trackno  = -1;
 427     }
 428   else if (ngtrack >= 1)
 429     {
 430       // one or more than one photon track + charged track
 431       // find out which track deposits maximum energy and
 432       // assign that track number and track pid
 433
 434       trenergy  = new Float_t [nghadtrack];
 435       trpid     = new Int_t [nghadtrack];
 436       sortcoord = new Int_t [2*nghadtrack];
 437       for (Int_t i = 0; i < ngtrack; i++)
 438         {
 439           trenergy[i] = 0.;
 440           trpid[i]    = -1;
 441           for (Int_t j = 0; j < nentry; j++)
 442             {
 443               if (ipid[j] == 22 && itra[j] == itra[phentry[i]])
 444                 {
 445                   trenergy[i] += cadc[j];
 446                   trpid[i]     = 22;
 447                 }
 448             }
 449         }
 450       for (Int_t i = ngtrack; i < nghadtrack; i++)
 451         {
 452           trenergy[i] = 0.;
 453           trpid[i]    = -1;
 454           for (Int_t j = 0; j < nentry; j++)
 455             {
 456               if (ipid[j] != 22 && itra[j] == itra[hadentry[i-ngtrack]])
 457                 {
 458                   trenergy[i] += cadc[j];
 459                   trpid[i]     = ipid[j];
 460                 }
 461             }
 462         }
 463
 464       Bool_t jsort = true;
 465       TMath::Sort(nghadtrack,trenergy,sortcoord,jsort);
 466
 467       Int_t gtr = sortcoord[0];
 468       if (trpid[gtr] == 22)
 469         {
 470           trackpid = 22;
 471           trackno  = itra[phentry[gtr]];   // highest adc track
 472         }
 473       else
 474         {
 475           trackpid = 8;
 476           trackno = -1;
 477         }
 478
 479       delete [] trenergy;
 480       delete [] trpid;
 481       delete [] sortcoord;
 482
 483     }   // end of ngtrack >= 1
 484
 485   delete [] phentry;
 486   delete [] hadentry;
 487
 488 }
 489 //--------------------------------------------------------------------//
 490 void AliPMDtracker::SetVertex(Double_t vtx[3], Double_t evtx[3])
 491 {
 492   fXvertex = vtx[0];
 493   fYvertex = vtx[1];
 494   fZvertex = vtx[2];
 495   fSigmaX  = evtx[0];
 496   fSigmaY  = evtx[1];
 497   fSigmaZ  = evtx[2];
 498 }
 499 //--------------------------------------------------------------------//
 500 void AliPMDtracker::ResetClusters()
 501 {
 502   if (fRecpoints) fRecpoints->Clear();
 503 }
 504 //--------------------------------------------------------------------//