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