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