PHOS/AliPHOSTrackSegmentMakerv1.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.92  2007/08/28 12:55:08  policheh
  21  * Loaders removed from the reconstruction code (C.Cheshkov)
  22  *
  23  * Revision 1.91  2007/08/07 14:12:03  kharlov
  24  * Quality assurance added (Yves Schutz)
  25  *
  26  * Revision 1.90  2007/07/11 13:43:30  hristov
  27  * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
  28  *
  29  * Revision 1.89  2007/07/03 08:13:04  kharlov
  30  * Bug fix in CPV local coordinates
  31  *
  32  * Revision 1.88  2007/06/27 09:11:07  kharlov
  33  * Bug fix for CPV-EMC distance
  34  *
  35  * Revision 1.87  2007/05/04 14:49:29  policheh
  36  * AliPHOSRecPoint inheritance from AliCluster
  37  *
  38  * Revision 1.86  2007/04/02 15:00:16  cvetan
  39  * No more calls to gAlice in the reconstruction
  40  *
  41  * Revision 1.85  2007/03/28 19:18:15  kharlov
  42  * RecPoints recalculation in TSM removed
  43  *
  44  * Revision 1.84  2007/03/07 07:01:21  hristov
  45  * Fixing copy/paste erro. Additional protections
  46  *
  47  * Revision 1.83  2007/03/06 21:07:37  kharlov
  48  * DP: xz CPV-EMC distance filled to TS
  49  *
  50  * Revision 1.82  2007/03/06 06:54:48  kharlov
  51  * DP:Calculation of cluster properties dep. on vertex added
  52  *
  53  * Revision 1.81  2007/02/05 10:02:40  kharlov
  54  * Module numbering is corrected
  55  *
  56  * Revision 1.80  2006/08/28 10:01:56  kharlov
  57  * Effective C++ warnings fixed (Timur Pocheptsov)
  58  *
  59  * Revision 1.79  2006/04/25 12:41:15  hristov
  60  * Moving non-persistent data to AliESDfriend (Yu.Belikov)
  61  *
  62  * Revision 1.78  2005/11/18 13:04:51  hristov
  63  * Bug fix
  64  *
  65  * Revision 1.77  2005/11/17 23:34:36  hristov
  66  * Corrected logics
  67  *
  68  * Revision 1.76  2005/11/17 22:29:12  hristov
  69  * Faster version, no attempt to match tracks outside the PHOS acceptance
  70  *
  71  * Revision 1.75  2005/11/17 12:35:27  hristov
  72  * Use references instead of objects. Avoid to create objects when they are not really needed
  73  *
  74  * Revision 1.74  2005/07/08 14:01:36  hristov
  75  * Tracking in non-uniform nmagnetic field (Yu.Belikov)
  76  *
  77  * Revision 1.73  2005/05/28 14:19:05  schutz
  78  * Compilation warnings fixed by T.P.
  79  *
  80  */
  81
  82 //_________________________________________________________________________
  83 // Implementation version 1 of algorithm class to construct PHOS track segments
  84 // Track segment for PHOS is list of
  85 //        EMC RecPoint + (possibly) CPV RecPoint
  86 // To find TrackSegments we do the following:
  87 //  for each EMC RecPoints we look at
  88 //   CPV RecPoints in the radious fRcpv.
  89 //  If there is such a CPV RecPoint,
  90 //   we make "Link" it is just indexes of EMC and CPV RecPoint and distance
  91 //   between them in the PHOS plane.
  92 //  Then we sort "Links" and starting from the
  93 //   least "Link" pointing to the unassined EMC and CPV RecPoints assing them to
  94 //   new TrackSegment.
  95 // If there is no CPV RecPoint we make TrackSegment
  96 // consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
  97 //// In principle this class should be called from AliPHOSReconstructor, but
  98 // one can use it as well in standalone mode.
  99 // Use  case:
 100 //  root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
 101 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
 102 //               // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
 103 //               // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")
 104 //  root [1] t->ExecuteTask()
 105 //  root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
 106 //  root [4] t->ExecuteTask("deb all time")
 107 //
 108 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
 109 //
 110
 111 // --- ROOT system ---
 112 #include "TTree.h"
 113 #include "TBenchmark.h"
 114
 115 // --- Standard library ---
 116 #include "Riostream.h"
 117 // --- AliRoot header files ---
 118 #include "AliPHOSGeometry.h"
 119 #include "AliPHOSTrackSegmentMakerv1.h"
 120 #include "AliPHOSTrackSegment.h"
 121 #include "AliPHOSLink.h"
 122 #include "AliESDEvent.h"
 123 #include "AliESDtrack.h"
 124 #include "AliPHOSQADataMaker.h"
 125 #include "AliPHOSEmcRecPoint.h"
 126 #include "AliPHOSCpvRecPoint.h"
 127
 128 ClassImp( AliPHOSTrackSegmentMakerv1)
 129
 130
 131 //____________________________________________________________________________
 132 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() :
 133   AliPHOSTrackSegmentMaker(),
 134   fDefaultInit(kTRUE),
 135   fWrite(kFALSE),
 136   fNTrackSegments(0),
 137   fRcpv(0.f),
 138   fRtpc(0.f),
 139   fVtx(0.f),
 140   fLinkUpArray(0),
 141   fEmcFirst(0),
 142   fEmcLast(0),
 143   fCpvFirst(0),
 144   fCpvLast(0),
 145   fModule(0),
 146   fTrackSegments(NULL)
 147 {
 148   // default ctor (to be used mainly by Streamer)
 149   InitParameters() ;
 150 }
 151
 152 //____________________________________________________________________________
 153 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(AliPHOSGeometry *geom) :
 154   AliPHOSTrackSegmentMaker(geom),
 155   fDefaultInit(kFALSE),
 156   fWrite(kFALSE),
 157   fNTrackSegments(0),
 158   fRcpv(0.f),
 159   fRtpc(0.f),
 160   fVtx(0.f),
 161   fLinkUpArray(0),
 162   fEmcFirst(0),
 163   fEmcLast(0),
 164   fCpvFirst(0),
 165   fCpvLast(0),
 166   fModule(0),
 167   fTrackSegments(NULL)
 168 {
 169   // ctor
 170   InitParameters() ;
 171   Init() ;
 172   fESD = 0;
 173 }
 174
 175
 176 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const AliPHOSTrackSegmentMakerv1 & tsm) :
 177   AliPHOSTrackSegmentMaker(tsm),
 178   fDefaultInit(kFALSE),
 179   fWrite(kFALSE),
 180   fNTrackSegments(0),
 181   fRcpv(0.f),
 182   fRtpc(0.f),
 183   fVtx(0.f),
 184   fLinkUpArray(0),
 185   fEmcFirst(0),
 186   fEmcLast(0),
 187   fCpvFirst(0),
 188   fCpvLast(0),
 189   fModule(0),
 190   fTrackSegments(NULL)
 191 {
 192   // cpy ctor: no implementation yet
 193   // requested by the Coding Convention
 194   Fatal("cpy ctor", "not implemented") ;
 195 }
 196
 197
 198 //____________________________________________________________________________
 199  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
 200 {
 201   // dtor
 202   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
 203   if (!fDefaultInit)
 204     delete fLinkUpArray ;
 205   if (fTrackSegments) {
 206     fTrackSegments->Delete();
 207     delete fTrackSegments;
 208   }
 209 }
 210
 211 //____________________________________________________________________________
 212 void  AliPHOSTrackSegmentMakerv1::FillOneModule()
 213 {
 214   // Finds first and last indexes between which
 215   // clusters from one PHOS module are
 216
 217   //First EMC clusters
 218   Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
 219   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
 220         ((dynamic_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
 221       fEmcLast ++)  ;
 222
 223   //Now CPV clusters
 224   Int_t totalCpv = fCPVRecPoints->GetEntriesFast() ;
 225
 226     for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) &&
 227          ((dynamic_cast<AliPHOSRecPoint *>(fCPVRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule );
 228        fCpvLast ++) ;
 229
 230 }
 231
 232 //____________________________________________________________________________
 233 void  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
 234                                                          AliPHOSCpvRecPoint * cpvClu,
 235                                                          Int_t &trackindex,
 236                                                          Float_t &dx, Float_t &dz) const
 237 {
 238   // Calculates the distance between the EMC RecPoint and the CPV RecPoint
 239   // Clusters are sorted in "rows" and "columns" of width 1 cm
 240
 241 //  Float_t delta = 1 ;  // Width of the rows in sorting of RecPoints (in cm)
 242 //                       // if you change this value, change it as well in xxxRecPoint::Compare()
 243   Float_t distance2Track = fRtpc ;
 244
 245   trackindex = -1 ; // closest track within fRCpv
 246
 247   TVector3 vecEmc ;   // Local position of EMC recpoint
 248   TVector3 vecP ;     // Momentum direction at CPV plain
 249   TVector3 vecPloc ;     // Momentum direction at CPV plain
 250
 251   //toofar = kTRUE ;
 252   if(emcClu->GetPHOSMod() != cpvClu->GetPHOSMod()){
 253     dx=999. ;
 254     dz=999. ;
 255     return ;
 256   }
 257
 258   emcClu->GetLocalPosition(vecEmc) ;
 259
 260   Double_t xCPV,zCPV ; //EMC-projected coordinates of CPV cluster
 261   TVector3 cpvGlobal; // Global position of the CPV recpoint
 262   fGeom->GetGlobalPHOS((AliPHOSRecPoint*)cpvClu,cpvGlobal);
 263   Double_t vtxCPV[3]={cpvGlobal.X(),cpvGlobal.Y(),cpvGlobal.Z()} ;
 264   Int_t dummyMod ;
 265
 266   if (fESD == 0x0) {
 267      //if no track information available, assume straight line from IP to emcal
 268      fGeom->ImpactOnEmc(vtxCPV,cpvGlobal.Theta(),cpvGlobal.Phi(),dummyMod,zCPV,xCPV) ;
 269      dx=xCPV - vecEmc.X() ;
 270      dz=zCPV - vecEmc.Z() ;
 271      return ;
 272   }
 273
 274   //if there is ESD try to correct distance using TPC information on particle direct in CPV
 275   if (fESD != 0x0) {
 276
 277     Double_t rCPV = cpvGlobal.Pt() ;// Radius from IP to current point
 278
 279     // Extrapolate the global track direction if any to CPV and find the closest track
 280     Int_t nTracks = fESD->GetNumberOfTracks();
 281     Int_t iClosestTrack = -1;
 282     Double_t minDistance = 1.e6;
 283     TVector3 inPHOS ;
 284
 285     AliESDtrack *track;
 286     Double_t xyz[3] ;
 287     Double_t pxyz[3];
 288     for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
 289       track = fESD->GetTrack(iTrack);
 290       if (!track->GetXYZAt(rCPV, fESD->GetMagneticField(), xyz))
 291            continue; //track coord on the cylinder of PHOS radius
 292       if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
 293            continue;
 294       //Check if this track hits PHOS
 295       inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
 296       distance2Track = inPHOS.Angle(cpvGlobal) ;
 297       // Find the closest track to the CPV recpoint
 298       if (distance2Track < minDistance) {
 299          minDistance = distance2Track;
 300          iClosestTrack = iTrack;
 301       }
 302     }
 303
 304     if (iClosestTrack != -1) {
 305       track = fESD->GetTrack(iClosestTrack);
 306       if (track->GetPxPyPzAt(rCPV, fESD->GetMagneticField(), pxyz)) { // track momentum ibid.
 307         vecP.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
 308         Int_t dummyMod ;
 309         fGeom->ImpactOnEmc(vtxCPV,vecP.Theta(),vecP.Phi(),dummyMod,zCPV,xCPV) ;
 310       }
 311     }
 312
 313     if(minDistance < fRtpc ){
 314       trackindex = iClosestTrack ;
 315     }
 316   }
 317   if(trackindex!=-1){
 318     // If the closest global track is found, calculate EMC-CPV distance from it
 319     dx=xCPV - vecEmc.X() ;
 320     dz=zCPV - vecEmc.Z() ;
 321   }
 322   else{
 323     // If no global track was found, just take the nearest CPV point
 324     fGeom->ImpactOnEmc(vtxCPV,cpvGlobal.Theta(),cpvGlobal.Phi(),dummyMod,zCPV,xCPV) ;
 325     dx=xCPV - vecEmc.X() ;
 326     dz=zCPV - vecEmc.Z() ;
 327   }
 328   return ;
 329 }
 330 //____________________________________________________________________________
 331 void  AliPHOSTrackSegmentMakerv1::Init()
 332 {
 333   // Make all memory allocations that are not possible in default constructor
 334
 335   fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000);
 336   fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100);
 337   fTrackSegments->SetName("TRACKS");
 338 }
 339
 340 //____________________________________________________________________________
 341 void  AliPHOSTrackSegmentMakerv1::InitParameters()
 342 {
 343   //Initializes parameters
 344   fRcpv      = 10. ;
 345   fRtpc      = 4. ;
 346   fEmcFirst  = 0 ;
 347   fEmcLast   = 0 ;
 348   fCpvFirst  = 0 ;
 349   fCpvLast   = 0 ;
 350   fLinkUpArray = 0 ;
 351   fWrite                   = kTRUE ;
 352 }
 353
 354
 355 //____________________________________________________________________________
 356 void  AliPHOSTrackSegmentMakerv1::MakeLinks()const
 357 {
 358   // Finds distances (links) between all EMC and CPV clusters,
 359   // which are not further apart from each other than fRcpv
 360   // and sort them in accordance with this distance
 361
 362   fLinkUpArray->Clear() ;
 363
 364   AliPHOSCpvRecPoint * cpv ;
 365   AliPHOSEmcRecPoint * emcclu ;
 366
 367   Int_t iLinkUp  = 0 ;
 368
 369   Int_t iEmcRP;
 370   for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
 371     emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP)) ;
 372
 373     //Bool_t toofar ;
 374     Int_t iCpv = 0 ;
 375     for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) {
 376
 377       cpv = dynamic_cast<AliPHOSCpvRecPoint *>(fCPVRecPoints->At(iCpv)) ;
 378       Int_t track = -1 ;
 379       Float_t dx,dz ;
 380       GetDistanceInPHOSPlane(emcclu, cpv, track,dx,dz) ;
 381       if(TMath::Sqrt(dx*dx+dz*dz) < fRcpv ){
 382         new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(dx, dz, iEmcRP, iCpv, track) ;
 383       }
 384     }
 385   }
 386
 387   fLinkUpArray->Sort() ;  //first links with smallest distances
 388 }
 389
 390 //____________________________________________________________________________
 391 void  AliPHOSTrackSegmentMakerv1::MakePairs()
 392 {
 393   // Using the previously made list of "links", we found the smallest link - i.e.
 394   // link with the least distance between EMC and CPV and pointing to still
 395   // unassigned RecParticles. We assign these RecPoints to TrackSegment and
 396   // remove them from the list of "unassigned".
 397
 398   //Make arrays to mark clusters already chosen
 399   Int_t * emcExist = 0;
 400   if(fEmcLast > fEmcFirst)
 401     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
 402
 403   Int_t index;
 404   for(index = 0; index <fEmcLast-fEmcFirst; index ++)
 405     emcExist[index] = 1 ;
 406
 407   Bool_t * cpvExist = 0;
 408   if(fCpvLast > fCpvFirst)
 409     cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
 410   for(index = 0; index <fCpvLast-fCpvFirst; index ++)
 411     cpvExist[index] = kTRUE ;
 412
 413
 414   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
 415   TIter nextUp(fLinkUpArray) ;
 416
 417   AliPHOSLink * linkUp ;
 418
 419   AliPHOSCpvRecPoint * nullpointer = 0 ;
 420
 421   while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){
 422
 423     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
 424
 425       if(cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist
 426          Float_t dx,dz ;
 427          linkUp->GetXZ(dx,dz) ;
 428          new ((* fTrackSegments)[fNTrackSegments])
 429            AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) ,
 430                                dynamic_cast<AliPHOSCpvRecPoint *>(fCPVRecPoints->At(linkUp->GetCpv())) ,
 431                                linkUp->GetTrack(),dx,dz) ;
 432
 433        (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 434        fNTrackSegments++ ;
 435        emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found
 436        //mark CPV recpoint as already used
 437        cpvExist[linkUp->GetCpv()-fCpvFirst] = kFALSE ;
 438       } //if CpvUp still exist
 439     }
 440   }
 441
 442   //look through emc recPoints left without CPV
 443   if(emcExist){ //if there is emc rec point
 444     Int_t iEmcRP ;
 445     for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
 446       if(emcExist[iEmcRP] > 0 ){
 447        new ((*fTrackSegments)[fNTrackSegments])
 448          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP+fEmcFirst)),
 449                            nullpointer) ;
 450        (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
 451        fNTrackSegments++;
 452       }
 453     }
 454   }
 455   delete [] emcExist ;
 456   delete [] cpvExist ;
 457 }
 458
 459 //____________________________________________________________________________
 460 void AliPHOSTrackSegmentMakerv1::Clusters2TrackSegments(Option_t *option)
 461 {
 462   // Steering method to perform track segment construction for the current event
 463   // Returns an array with the found track-segments.
 464
 465   if(strstr(option,"tim"))
 466     gBenchmark->Start("PHOSTSMaker");
 467
 468   if(strstr(option,"print")) {
 469     Print() ;
 470     return ;
 471   }
 472
 473   //Make some initializations
 474   fNTrackSegments = 0 ;
 475   fEmcFirst = 0 ;
 476   fEmcLast  = 0 ;
 477   fCpvFirst = 0 ;
 478   fCpvLast  = 0 ;
 479
 480   fTrackSegments->Clear();
 481
 482   //   if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
 483
 484   for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) {
 485     FillOneModule() ;
 486     MakeLinks() ;
 487     MakePairs() ;
 488   }
 489
 490   if(strstr(option,"deb"))
 491     PrintTrackSegments(option);
 492
 493   if(strstr(option,"tim")){
 494     gBenchmark->Stop("PHOSTSMaker");
 495     Info("Exec", "took %f seconds for making TS",
 496          gBenchmark->GetCpuTime("PHOSTSMaker"));
 497   }
 498 }
 499
 500 //____________________________________________________________________________
 501 void AliPHOSTrackSegmentMakerv1::Print(const Option_t *)const
 502 {
 503   //  Print TrackSegmentMaker parameters
 504
 505   TString message("") ;
 506   if( strcmp(GetName(), "") != 0 ) {
 507     message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ;
 508     message += "Making Track segments\n" ;
 509     message += "with parameters:\n" ;
 510     message += "     Maximal EMC - CPV distance (cm) %f\n" ;
 511     message += "============================================\n" ;
 512     Info("Print", message.Data(),fRcpv) ;
 513   }
 514   else
 515     Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
 516 }
 517
 518 //____________________________________________________________________________
 519 void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
 520 {
 521   // option deb - prints # of found TrackSegments
 522   // option deb all - prints as well indexed of found RecParticles assigned to the TS
 523
 524   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
 525   printf("        Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() );
 526
 527   if(strstr(option,"all")) {  // printing found TS
 528     printf("TrackSegment #  EMC RP#  CPV RP#\n") ;
 529     Int_t index;
 530     for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
 531       AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ;
 532       printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ;
 533     }
 534   }
 535 }