PHOS/AliPHOSTrackSegmentMakerv2.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 /* $Id$ */
  16
  17 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.6  2007/08/22 09:20:50  hristov
  21  * Updated QA classes (Yves)
  22  *
  23  * Revision 1.5  2007/07/11 13:43:30  hristov
  24  * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
  25  *
  26  * Revision 1.4  2007/05/04 14:49:29  policheh
  27  * AliPHOSRecPoint inheritance from AliCluster
  28  *
  29  * Revision 1.3  2007/04/25 19:39:42  kharlov
  30  * Track extracpolation improved
  31  *
  32  * Revision 1.2  2007/04/01 19:16:52  kharlov
  33  * D.P.: Produce EMCTrackSegments using TPC/ITS tracks (no CPV)
  34  *
  35  *
  36  */
  37
  38 //_________________________________________________________________________
  39 // Implementation version 2 of algorithm class to construct PHOS track segments
  40 // Track segment for PHOS is list of
  41 //        EMC RecPoint + (possibly) projection of TPC track
  42 // To find TrackSegments we do the following:
  43 //  for each EMC RecPoints we look at
  44 //   TPC projections radius fRtpc.
  45 //  If there is such a track
  46 //   we make "Link" it is just indexes of EMC and TPC track and distance
  47 //   between them in the PHOS plane.
  48 //  Then we sort "Links" and starting from the
  49 //   least "Link" pointing to the unassined EMC and TPC assing them to
  50 //   new TrackSegment.
  51 // If there is no TPC track we make TrackSegment
  52 // consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
  53 //// In principle this class should be called from AliPHOSReconstructor, but
  54 // one can use it as well in standalone mode.
  55 // Use  case:
  56 //  root [0] AliPHOSTrackSegmentMakerv2 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
  57 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  58 //               // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
  59 //               // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")
  60 //  root [1] t->ExecuteTask()
  61 //  root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
  62 //  root [4] t->ExecuteTask("deb all time")
  63 //
  64 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
  65 //
  66
  67 // --- ROOT system ---
  68 #include "TFile.h"
  69 #include "TTree.h"
  70 #include "TBenchmark.h"
  71
  72 // --- Standard library ---
  73 #include "Riostream.h"
  74 // --- AliRoot header files ---
  75 #include "AliPHOSGeometry.h"
  76 #include "AliPHOSTrackSegmentMakerv2.h"
  77 #include "AliPHOSTrackSegment.h"
  78 #include "AliPHOSLink.h"
  79 #include "AliPHOSGetter.h"
  80
  81 #include "AliESDEvent.h"
  82 #include "AliESDtrack.h"
  83
  84 ClassImp( AliPHOSTrackSegmentMakerv2)
  85
  86
  87 //____________________________________________________________________________
  88 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2() :
  89   AliPHOSTrackSegmentMaker(),
  90   fDefaultInit(kTRUE),
  91   fWrite(kFALSE),
  92   fNTrackSegments(0),
  93   fRtpc(0.f),
  94   fVtx(0.,0.,0.),
  95   fLinkUpArray(0),
  96   fTPCtracks(),
  97   fEmcFirst(0),
  98   fEmcLast(0),
  99   fModule(0),
 100   fTrackSegments(NULL)
 101 {
 102   // default ctor (to be used mainly by Streamer)
 103   InitParameters() ;
 104 }
 105
 106 //____________________________________________________________________________
 107 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(AliPHOSGeometry *geom) :
 108   AliPHOSTrackSegmentMaker(geom),
 109   fDefaultInit(kFALSE),
 110   fWrite(kFALSE),
 111   fNTrackSegments(0),
 112   fRtpc(0.f),
 113   fVtx(0.,0.,0.),
 114   fLinkUpArray(0),
 115   fTPCtracks(),
 116   fEmcFirst(0),
 117   fEmcLast(0),
 118   fModule(0),
 119   fTrackSegments(NULL)
 120 {
 121   // ctor
 122   InitParameters() ;
 123   Init() ;
 124   fESD = 0;
 125 }
 126
 127 //____________________________________________________________________________
 128 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
 129   AliPHOSTrackSegmentMaker(tsm),
 130   fDefaultInit(kFALSE),
 131   fWrite(kFALSE),
 132   fNTrackSegments(0),
 133   fRtpc(0.f),
 134   fVtx(0.,0.,0.),
 135   fLinkUpArray(0),
 136   fTPCtracks(),
 137   fEmcFirst(0),
 138   fEmcLast(0),
 139   fModule(0),
 140   fTrackSegments(NULL)
 141 {
 142   // cpy ctor: no implementation yet
 143   // requested by the Coding Convention
 144   Fatal("cpy ctor", "not implemented") ;
 145 }
 146
 147
 148 //____________________________________________________________________________
 149  AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
 150 {
 151   // dtor
 152   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
 153   if (!fDefaultInit)
 154     delete fLinkUpArray ;
 155   if (fTrackSegments) {
 156     fTrackSegments->Delete();
 157     delete fTrackSegments;
 158   }
 159 }
 160
 161 //____________________________________________________________________________
 162 void  AliPHOSTrackSegmentMakerv2::FillOneModule()
 163 {
 164   // Finds first and last indexes between which
 165   // clusters from one PHOS module are
 166
 167   //First EMC clusters
 168   Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
 169   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
 170         ((dynamic_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
 171       fEmcLast ++)  ;
 172
 173   //Now TPC tracks
 174   if(fESD){
 175     //Do it ones, only first time
 176     if(fModule==1){
 177       Int_t nTracks = fESD->GetNumberOfTracks();
 178
 179       Int_t nPHOSmod = fGeom->GetNModules() ;
 180       if(fTPCtracks[0].size()<(UInt_t)nTracks){
 181         for(Int_t i=0; i<nPHOSmod; i++)
 182           fTPCtracks[i].resize(nTracks) ;
 183       }
 184       for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
 185       TVector3 inPHOS ;
 186
 187       //In this particular case we use fixed vertex position at zero
 188       Double_t vtx[3]={0.,0.,0.} ;
 189       AliESDtrack *track;
 190       Double_t xyz[3] ;
 191       Double_t rEMC = fGeom->GetIPtoCrystalSurface() ; //Use here ideal geometry
 192       for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
 193         track = fESD->GetTrack(iTrack);
 194         for(Int_t iTestMod=1; iTestMod<=nPHOSmod; iTestMod++){
 195            Double_t modAngle=270.+fGeom->GetPHOSAngle(iTestMod) ;
 196            modAngle*=TMath::Pi()/180. ;
 197            track->Rotate(modAngle);
 198            if (!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz))
 199              continue; //track coord on the cylinder of PHOS radius
 200            if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
 201              continue;
 202            //Check if this track hits PHOS
 203            inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
 204            Int_t modNum ;
 205            Double_t x,z ;
 206            fGeom->ImpactOnEmc(vtx, inPHOS.Theta(), inPHOS.Phi(), modNum, z, x) ;
 207            if(modNum==iTestMod){
 208              //Mark this track as one belonging to module
 209              TrackInPHOS_t &t = fTPCtracks[modNum-1][fNtpcTracks[modNum-1]++] ;
 210              t.track = track ;
 211              t.x = x ;
 212              t.z = z ;
 213              break ;
 214            }
 215         }
 216       }
 217     }
 218   }
 219
 220 }
 221
 222 //____________________________________________________________________________
 223 void  AliPHOSTrackSegmentMakerv2::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
 224                                                          AliESDtrack *track,
 225                                                          Float_t &dx, Float_t &dz) const
 226 {
 227   // Calculates the distance between the EMC RecPoint and the CPV RecPoint
 228   // Clusters are sorted in "rows" and "columns" of width 1 cm
 229
 230   TVector3 emcGlobal; // Global position of the EMC recpoint
 231   fGeom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);
 232
 233   //Calculate actual distance to the PHOS surface
 234   Double_t modAngle=270.+fGeom->GetPHOSAngle(emcClu->GetPHOSMod()) ;
 235   modAngle*=TMath::Pi()/180. ;
 236   Double_t rEMC = emcGlobal.X()*TMath::Cos(modAngle)+emcGlobal.Y()*TMath::Sin(modAngle) ;
 237   track->Rotate(modAngle);
 238   Double_t xyz[3] ;
 239   if(!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)){
 240     dx=999. ;
 241     dz=999. ;
 242     return ; //track coord on the cylinder of PHOS radius
 243   }
 244   if((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0){
 245     dx=999. ;
 246     dz=999. ;
 247     return ;
 248   }
 249
 250   dx=(emcGlobal.X()-xyz[0])*TMath::Sin(modAngle)-(emcGlobal.Y()-xyz[1])*TMath::Cos(modAngle) ;
 251   dx=TMath::Sign(dx,(Float_t)(emcGlobal.X()-xyz[0])) ; //set direction
 252   dz=emcGlobal.Z()-xyz[2] ;
 253
 254   return ;
 255 }
 256 //____________________________________________________________________________
 257 void  AliPHOSTrackSegmentMakerv2::Init()
 258 {
 259   // Make all memory allocations that are not possible in default constructor
 260
 261   fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000);
 262   fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100);
 263 }
 264
 265 //____________________________________________________________________________
 266 void  AliPHOSTrackSegmentMakerv2::InitParameters()
 267 {
 268   //Initializes parameters
 269   fRtpc      = 4. ;
 270   fEmcFirst  = 0 ;
 271   fEmcLast   = 0 ;
 272   fLinkUpArray = 0 ;
 273   fWrite                   = kTRUE ;
 274 }
 275
 276
 277 //____________________________________________________________________________
 278 void  AliPHOSTrackSegmentMakerv2::MakeLinks()
 279 {
 280   // Finds distances (links) between all EMC and CPV clusters,
 281   // which are not further apart from each other than fRcpv
 282   // and sort them in accordance with this distance
 283
 284   fLinkUpArray->Clear() ;
 285
 286   AliPHOSEmcRecPoint * emcclu ;
 287
 288   Int_t iLinkUp  = 0 ;
 289
 290   Int_t iEmcRP;
 291   for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
 292     emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP)) ;
 293     TVector3 vecEmc ;
 294     emcclu->GetLocalPosition(vecEmc) ;
 295     Int_t mod=emcclu->GetPHOSMod() ;
 296     for(Int_t itr=0; itr<fNtpcTracks[mod-1] ; itr++){
 297       TrackInPHOS_t &t = fTPCtracks[mod-1][itr] ;
 298       //calculate raw distance
 299       Double_t rawdx = t.x - vecEmc.X() ;
 300       Double_t rawdz = t.z - vecEmc.Z() ;
 301       if(TMath::Sqrt(rawdx*rawdx+rawdz*rawdz)<fRtpc){
 302         Float_t dx,dz ;
 303         //calculate presize difference accounting misalignment
 304         GetDistanceInPHOSPlane(emcclu, t.track, dx,dz) ;
 305         Int_t itrack = t.track->GetID() ;
 306         new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(dx, dz, iEmcRP, itrack, -1) ;
 307       }
 308     }
 309   }
 310   fLinkUpArray->Sort() ;  //first links with smallest distances
 311 }
 312
 313 //____________________________________________________________________________
 314 void  AliPHOSTrackSegmentMakerv2::MakePairs()
 315 {
 316   // Using the previously made list of "links", we found the smallest link - i.e.
 317   // link with the least distance between EMC and CPV and pointing to still
 318   // unassigned RecParticles. We assign these RecPoints to TrackSegment and
 319   // remove them from the list of "unassigned".
 320
 321   //Make arrays to mark clusters already chosen
 322   Int_t * emcExist = 0;
 323   if(fEmcLast > fEmcFirst)
 324     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
 325
 326   Int_t index;
 327   for(index = 0; index <fEmcLast-fEmcFirst; index ++)
 328     emcExist[index] = 1 ;
 329
 330
 331   Bool_t * tpcExist = 0;
 332   Int_t nTracks = fESD->GetNumberOfTracks();
 333   if(nTracks>0)
 334     tpcExist = new Bool_t[nTracks] ;
 335
 336   for(index = 0; index <nTracks; index ++)
 337     tpcExist[index] = 1 ;
 338
 339
 340   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
 341   TIter nextUp(fLinkUpArray) ;
 342
 343   AliPHOSLink * linkUp ;
 344
 345   AliPHOSCpvRecPoint * nullpointer = 0 ;
 346
 347   while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){
 348
 349     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
 350
 351       if(tpcExist[linkUp->GetCpv()]){ //Track still exist
 352          Float_t dx,dz ;
 353          linkUp->GetXZ(dx,dz) ;
 354          new ((* fTrackSegments)[fNTrackSegments])
 355            AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) ,
 356                                nullpointer,
 357                                linkUp->GetTrack(),dx,dz) ;
 358        (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 359        fNTrackSegments++ ;
 360        emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found
 361        //mark track as already used
 362        tpcExist[linkUp->GetCpv()] = kFALSE ;
 363       } //if CpvUp still exist
 364     }
 365   }
 366
 367   //look through emc recPoints left without CPV
 368   if(emcExist){ //if there is emc rec point
 369     Int_t iEmcRP ;
 370     for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
 371       if(emcExist[iEmcRP] > 0 ){
 372        new ((*fTrackSegments)[fNTrackSegments])
 373          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP+fEmcFirst)),
 374                            nullpointer) ;
 375        (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 376        fNTrackSegments++;
 377       }
 378     }
 379   }
 380   delete [] emcExist ;
 381   delete [] tpcExist ;
 382 }
 383
 384 //____________________________________________________________________________
 385 void  AliPHOSTrackSegmentMakerv2::Clusters2TrackSegments(Option_t *option)
 386 {
 387   // Steering method to perform track segment construction for events
 388   // in the range from fFirstEvent to fLastEvent.
 389   // This range is optionally set by SetEventRange().
 390   // if fLastEvent=-1 (by default), then process events until the end.
 391
 392   if(strstr(option,"tim"))
 393     gBenchmark->Start("PHOSTSMaker");
 394
 395   if(strstr(option,"print")) {
 396     Print() ;
 397     return ;
 398   }
 399
 400   //Make some initializations
 401   fNTrackSegments = 0 ;
 402   fEmcFirst = 0 ;
 403   fEmcLast  = 0 ;
 404
 405   fTrackSegments->Clear();
 406
 407   //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
 408
 409   for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) {
 410     FillOneModule() ;
 411     MakeLinks() ;
 412     MakePairs() ;
 413   }
 414
 415   if(strstr(option,"deb"))
 416     PrintTrackSegments(option);
 417
 418   if(strstr(option,"tim")){
 419     gBenchmark->Stop("PHOSTSMaker");
 420     Info("Exec", "took %f seconds for making TS",
 421          gBenchmark->GetCpuTime("PHOSTSMaker"));
 422   }
 423 }
 424
 425 //____________________________________________________________________________
 426 void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
 427 {
 428   //  Print TrackSegmentMaker parameters
 429
 430   TString message("") ;
 431   if( strcmp(GetName(), "") != 0 ) {
 432     message = "\n======== AliPHOSTrackSegmentMakerv2 ========\n" ;
 433     message += "Making Track segments\n" ;
 434     message += "with parameters:\n" ;
 435     message += "     Maximal EMC - TPC distance (cm) %f\n" ;
 436     message += "============================================\n" ;
 437     Info("Print", message.Data(),fRtpc) ;
 438   }
 439   else
 440     Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;
 441
 442 }
 443
 444 //____________________________________________________________________________
 445 void AliPHOSTrackSegmentMakerv2::PrintTrackSegments(Option_t * option)
 446 {
 447   // option deb - prints # of found TrackSegments
 448   // option deb all - prints as well indexed of found RecParticles assigned to the TS
 449
 450   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
 451   printf("        Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() );
 452
 453   if(strstr(option,"all")) {  // printing found TS
 454     printf("TrackSegment #  EMC RP#  CPV RP#\n") ;
 455     Int_t index;
 456     for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
 457       AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ;
 458       printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetTrackIndex() ) ;
 459     }
 460   }
 461 }