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