[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv2.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id$ */

/* History of cvs commits:
 *
 * $Log$
 * Revision 1.6  2007/08/22 09:20:50  hristov
 * Updated QA classes (Yves)
 *
 * Revision 1.5  2007/07/11 13:43:30  hristov
 * New class AliESDEvent, backward compatibility with the old AliESD (Christian)
 *
 * Revision 1.4  2007/05/04 14:49:29  policheh
 * AliPHOSRecPoint inheritance from AliCluster
 *
 * Revision 1.3  2007/04/25 19:39:42  kharlov
 * Track extracpolation improved
 *
 * Revision 1.2  2007/04/01 19:16:52  kharlov
 * D.P.: Produce EMCTrackSegments using TPC/ITS tracks (no CPV)
 *
 *
 */

//_________________________________________________________________________
// Implementation version 2 of algorithm class to construct PHOS track segments
// Track segment for PHOS is list of 
//        EMC RecPoint + (possibly) projection of TPC track
// To find TrackSegments we do the following: 
//  for each EMC RecPoints we look at
//   TPC projections radius fRtpc. 
//  If there is such a track
//   we make "Link" it is just indexes of EMC and TPC track and distance
//   between them in the PHOS plane. 
//  Then we sort "Links" and starting from the 
//   least "Link" pointing to the unassined EMC and TPC assing them to 
//   new TrackSegment. 
// If there is no TPC track we make TrackSegment 
// consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
//// In principle this class should be called from AliPHOSReconstructor, but 
// one can use it as well in standalone mode.
// Use  case:
//  root [0] AliPHOSTrackSegmentMakerv2 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
//  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
//               // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
//               // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")                       
//  root [1] t->ExecuteTask()
//  root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
//  root [4] t->ExecuteTask("deb all time") 
//                 
//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) 
//

// --- ROOT system ---
#include "TFile.h"
#include "TTree.h"
#include "TBenchmark.h"

// --- Standard library ---
#include "Riostream.h"
// --- AliRoot header files ---
#include "AliPHOSGeometry.h"
#include "AliPHOSTrackSegmentMakerv2.h"
#include "AliPHOSTrackSegment.h"
#include "AliPHOSLink.h"
#include "AliPHOSGetter.h"

#include "AliESDEvent.h"
#include "AliESDtrack.h"

ClassImp( AliPHOSTrackSegmentMakerv2) 


//____________________________________________________________________________
AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2() :
  AliPHOSTrackSegmentMaker(),
  fDefaultInit(kTRUE),
  fWrite(kFALSE),
  fNTrackSegments(0),
  fRtpc(0.f),
  fVtx(0.,0.,0.),
  fLinkUpArray(0),
  fTPCtracks(),
  fEmcFirst(0),
  fEmcLast(0),
  fModule(0),
  fTrackSegments(NULL)
{
  // default ctor (to be used mainly by Streamer)
  InitParameters() ; 
}

//____________________________________________________________________________
AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(AliPHOSGeometry *geom) :
  AliPHOSTrackSegmentMaker(geom),
  fDefaultInit(kFALSE),
  fWrite(kFALSE),
  fNTrackSegments(0),
  fRtpc(0.f),
  fVtx(0.,0.,0.),
  fLinkUpArray(0),
  fTPCtracks(),
  fEmcFirst(0),
  fEmcLast(0),
  fModule(0),
  fTrackSegments(NULL)
{
  // ctor
  InitParameters() ; 
  Init() ;
  fESD = 0;
}

//____________________________________________________________________________
AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
  AliPHOSTrackSegmentMaker(tsm),
  fDefaultInit(kFALSE),
  fWrite(kFALSE),
  fNTrackSegments(0),
  fRtpc(0.f),
  fVtx(0.,0.,0.),
  fLinkUpArray(0),
  fTPCtracks(),
  fEmcFirst(0),
  fEmcLast(0),
  fModule(0),
  fTrackSegments(NULL)
{
  // cpy ctor: no implementation yet
  // requested by the Coding Convention
  Fatal("cpy ctor", "not implemented") ;
}


//____________________________________________________________________________
 AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
{ 
  // dtor
  // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
  if (!fDefaultInit)  
    delete fLinkUpArray ;
  if (fTrackSegments) {
    fTrackSegments->Delete();
    delete fTrackSegments;
  }
}

//____________________________________________________________________________
void  AliPHOSTrackSegmentMakerv2::FillOneModule()
{
  // Finds first and last indexes between which 
  // clusters from one PHOS module are

  //First EMC clusters
  Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
  for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
	((dynamic_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
      fEmcLast ++)  ;
  
  //Now TPC tracks
  if(fESD){
    //Do it ones, only first time
    if(fModule==1){
      Int_t nTracks = fESD->GetNumberOfTracks();

      Int_t nPHOSmod = fGeom->GetNModules() ;
      if(fTPCtracks[0].size()<(UInt_t)nTracks){
        for(Int_t i=0; i<nPHOSmod; i++) 
          fTPCtracks[i].resize(nTracks) ;
      }
      for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
      TVector3 inPHOS ;
   
      //In this particular case we use fixed vertex position at zero
      Double_t vtx[3]={0.,0.,0.} ;
      AliESDtrack *track;
      Double_t xyz[3] ;
      Double_t rEMC = fGeom->GetIPtoCrystalSurface() ; //Use here ideal geometry 
      for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
        track = fESD->GetTrack(iTrack);
        for(Int_t iTestMod=1; iTestMod<=nPHOSmod; iTestMod++){
           Double_t modAngle=270.+fGeom->GetPHOSAngle(iTestMod) ;
           modAngle*=TMath::Pi()/180. ;
           track->Rotate(modAngle);  
           if (!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz))
             continue; //track coord on the cylinder of PHOS radius
           if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
             continue;
           //Check if this track hits PHOS
           inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
           Int_t modNum ; 
           Double_t x,z ;
           fGeom->ImpactOnEmc(vtx, inPHOS.Theta(), inPHOS.Phi(), modNum, z, x) ;
           if(modNum==iTestMod){
             //Mark this track as one belonging to module
             TrackInPHOS_t &t = fTPCtracks[modNum-1][fNtpcTracks[modNum-1]++] ; 
             t.track = track ;
             t.x = x ;
             t.z = z ;
             break ;
           }
        }
      }
    }
  } 

}

//____________________________________________________________________________
void  AliPHOSTrackSegmentMakerv2::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
                                                         AliESDtrack *track,
                                                         Float_t &dx, Float_t &dz) const
{
  // Calculates the distance between the EMC RecPoint and the CPV RecPoint
  // Clusters are sorted in "rows" and "columns" of width 1 cm

  TVector3 emcGlobal; // Global position of the EMC recpoint
  fGeom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);

  //Calculate actual distance to the PHOS surface
  Double_t modAngle=270.+fGeom->GetPHOSAngle(emcClu->GetPHOSMod()) ;
  modAngle*=TMath::Pi()/180. ;
  Double_t rEMC = emcGlobal.X()*TMath::Cos(modAngle)+emcGlobal.Y()*TMath::Sin(modAngle) ;
  track->Rotate(modAngle);
  Double_t xyz[3] ;
  if(!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)){
    dx=999. ;
    dz=999. ;
    return ; //track coord on the cylinder of PHOS radius
  }
  if((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0){
    dx=999. ; 
    dz=999. ;
    return ;
  } 

  dx=(emcGlobal.X()-xyz[0])*TMath::Sin(modAngle)-(emcGlobal.Y()-xyz[1])*TMath::Cos(modAngle) ;
  dx=TMath::Sign(dx,(Float_t)(emcGlobal.X()-xyz[0])) ; //set direction
  dz=emcGlobal.Z()-xyz[2] ;

  return ;
}
//____________________________________________________________________________
void  AliPHOSTrackSegmentMakerv2::Init()
{
  // Make all memory allocations that are not possible in default constructor
  
  fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
  fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100);
}

//____________________________________________________________________________
void  AliPHOSTrackSegmentMakerv2::InitParameters()
{
  //Initializes parameters
  fRtpc      = 4. ;
  fEmcFirst  = 0 ;    
  fEmcLast   = 0 ;   
  fLinkUpArray = 0 ;
  fWrite                   = kTRUE ;
}


//____________________________________________________________________________
void  AliPHOSTrackSegmentMakerv2::MakeLinks()
{ 
  // Finds distances (links) between all EMC and CPV clusters, 
  // which are not further apart from each other than fRcpv 
  // and sort them in accordance with this distance
  
  fLinkUpArray->Clear() ;    

  AliPHOSEmcRecPoint * emcclu ;

  Int_t iLinkUp  = 0 ;
  
  Int_t iEmcRP;
  for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
    emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP)) ;
    TVector3 vecEmc ;
    emcclu->GetLocalPosition(vecEmc) ;
    Int_t mod=emcclu->GetPHOSMod() ;
    for(Int_t itr=0; itr<fNtpcTracks[mod-1] ; itr++){
      TrackInPHOS_t &t = fTPCtracks[mod-1][itr] ;
      //calculate raw distance
      Double_t rawdx = t.x - vecEmc.X() ;
      Double_t rawdz = t.z - vecEmc.Z() ;
      if(TMath::Sqrt(rawdx*rawdx+rawdz*rawdz)<fRtpc){
        Float_t dx,dz ;
        //calculate presize difference accounting misalignment
        GetDistanceInPHOSPlane(emcclu, t.track, dx,dz) ;     
        Int_t itrack = t.track->GetID() ;
        new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(dx, dz, iEmcRP, itrack, -1) ;
      }      
    }
  } 
  fLinkUpArray->Sort() ;  //first links with smallest distances
}

//____________________________________________________________________________
void  AliPHOSTrackSegmentMakerv2::MakePairs()
{ 
  // Using the previously made list of "links", we found the smallest link - i.e. 
  // link with the least distance between EMC and CPV and pointing to still 
  // unassigned RecParticles. We assign these RecPoints to TrackSegment and 
  // remove them from the list of "unassigned". 

  //Make arrays to mark clusters already chosen
  Int_t * emcExist = 0;
  if(fEmcLast > fEmcFirst)
    emcExist = new Int_t[fEmcLast-fEmcFirst] ;
  
  Int_t index;
  for(index = 0; index <fEmcLast-fEmcFirst; index ++)
    emcExist[index] = 1 ;
  
  
  Bool_t * tpcExist = 0;
  Int_t nTracks = fESD->GetNumberOfTracks();
  if(nTracks>0)
    tpcExist = new Bool_t[nTracks] ;
  
  for(index = 0; index <nTracks; index ++)
    tpcExist[index] = 1 ;
  
  
  // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
  TIter nextUp(fLinkUpArray) ;
  
  AliPHOSLink * linkUp ;
  
  AliPHOSCpvRecPoint * nullpointer = 0 ;
  
  while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){  

    if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){

      if(tpcExist[linkUp->GetCpv()]){ //Track still exist
         Float_t dx,dz ; 
         linkUp->GetXZ(dx,dz) ;
	 new ((* fTrackSegments)[fNTrackSegments]) 
	   AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) , 
			       nullpointer,
			       linkUp->GetTrack(),dx,dz) ;
       (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
       fNTrackSegments++ ;
       emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found 
       //mark track as already used 
       tpcExist[linkUp->GetCpv()] = kFALSE ;
      } //if CpvUp still exist
    } 
  }        

  //look through emc recPoints left without CPV
  if(emcExist){ //if there is emc rec point
    Int_t iEmcRP ;
    for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst  ; iEmcRP++ ){
      if(emcExist[iEmcRP] > 0 ){
       new ((*fTrackSegments)[fNTrackSegments])  
         AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP+fEmcFirst)), 
                           nullpointer) ;
       (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
       fNTrackSegments++;    
      } 
    }
  }
  delete [] emcExist ; 
  delete [] tpcExist ; 
}

//____________________________________________________________________________
void  AliPHOSTrackSegmentMakerv2::Clusters2TrackSegments(Option_t *option)
{
  // Steering method to perform track segment construction for events
  // in the range from fFirstEvent to fLastEvent.
  // This range is optionally set by SetEventRange().
  // if fLastEvent=-1 (by default), then process events until the end.
  
  if(strstr(option,"tim"))
    gBenchmark->Start("PHOSTSMaker");
 
  if(strstr(option,"print")) {
    Print() ; 
    return ; 
  }
  
  //Make some initializations 
  fNTrackSegments = 0 ;
  fEmcFirst = 0 ;    
  fEmcLast  = 0 ;   
    
  fTrackSegments->Clear();

  //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
    
  for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) {
    FillOneModule() ; 
    MakeLinks() ;
    MakePairs() ;
  }

  if(strstr(option,"deb"))
    PrintTrackSegments(option);
    
  if(strstr(option,"tim")){
    gBenchmark->Stop("PHOSTSMaker");
    Info("Exec", "took %f seconds for making TS", 
	 gBenchmark->GetCpuTime("PHOSTSMaker")); 
  }
}

//____________________________________________________________________________
void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
{
  //  Print TrackSegmentMaker parameters

  TString message("") ;
  if( strcmp(GetName(), "") != 0 ) {
    message = "\n======== AliPHOSTrackSegmentMakerv2 ========\n" ; 
    message += "Making Track segments\n" ;
    message += "with parameters:\n" ; 
    message += "     Maximal EMC - TPC distance (cm) %f\n" ;
    message += "============================================\n" ;
    Info("Print", message.Data(),fRtpc) ;
  }
  else
    Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;

}

//____________________________________________________________________________
void AliPHOSTrackSegmentMakerv2::PrintTrackSegments(Option_t * option)
{
  // option deb - prints # of found TrackSegments
  // option deb all - prints as well indexed of found RecParticles assigned to the TS

  Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
  printf("        Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() ); 
  
  if(strstr(option,"all")) {  // printing found TS
    printf("TrackSegment #  EMC RP#  CPV RP#\n") ; 
    Int_t index;
    for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
      AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ; 
      printf("   %d           %d        %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetTrackIndex() ) ; 
    }	
  }
}
Commit	Line	Data
861b23c3	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$
9a2cdbdf	20	* Revision 1.6 2007/08/22 09:20:50 hristov
	21	* Updated QA classes (Yves)
	22	*
a5fa6165	23	* Revision 1.5 2007/07/11 13:43:30 hristov
	24	* New class AliESDEvent, backward compatibility with the old AliESD (Christian)
	25	*
af885e0f	26	* Revision 1.4 2007/05/04 14:49:29 policheh
	27	* AliPHOSRecPoint inheritance from AliCluster
	28	*
9ee9f78d	29	* Revision 1.3 2007/04/25 19:39:42 kharlov
	30	* Track extracpolation improved
	31	*
c978a69b	32	* Revision 1.2 2007/04/01 19:16:52 kharlov
	33	* D.P.: Produce EMCTrackSegments using TPC/ITS tracks (no CPV)
	34	*
861b23c3	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 ---
4b598bb1	68	#include "TFile.h"
861b23c3	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"
a5fa6165	80
af885e0f	81	#include "AliESDEvent.h"
861b23c3	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),
c978a69b	94	fVtx(0.,0.,0.),
861b23c3	95	fLinkUpArray(0),
c978a69b	96	fTPCtracks(),
861b23c3	97	fEmcFirst(0),
	98	fEmcLast(0),
	99	fModule(0),
9a2cdbdf	100	fTrackSegments(NULL)
861b23c3	101	{
	102	// default ctor (to be used mainly by Streamer)
	103	InitParameters() ;
861b23c3	104	}
	105
	106	//____________________________________________________________________________
9a2cdbdf	107	AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(AliPHOSGeometry *geom) :
9a2cdbdf	108	AliPHOSTrackSegmentMaker(geom),
861b23c3	109	fDefaultInit(kFALSE),
	110	fWrite(kFALSE),
	111	fNTrackSegments(0),
	112	fRtpc(0.f),
c978a69b	113	fVtx(0.,0.,0.),
861b23c3	114	fLinkUpArray(0),
c978a69b	115	fTPCtracks(),
861b23c3	116	fEmcFirst(0),
	117	fEmcLast(0),
	118	fModule(0),
9a2cdbdf	119	fTrackSegments(NULL)
861b23c3	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),
c978a69b	134	fVtx(0.,0.,0.),
861b23c3	135	fLinkUpArray(0),
c978a69b	136	fTPCtracks(),
861b23c3	137	fEmcFirst(0),
	138	fEmcLast(0),
	139	fModule(0),
9a2cdbdf	140	fTrackSegments(NULL)
861b23c3	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 ;
9a2cdbdf	155	if (fTrackSegments) {
	156	fTrackSegments->Delete();
	157	delete fTrackSegments;
	158	}
861b23c3	159	}
	160
	161	//____________________________________________________________________________
	162	void AliPHOSTrackSegmentMakerv2::FillOneModule()
	163	{
	164	// Finds first and last indexes between which
	165	// clusters from one PHOS module are
	166
861b23c3	167	//First EMC clusters
9a2cdbdf	168	Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
861b23c3	169	for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
9a2cdbdf	170	((dynamic_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
861b23c3	171	fEmcLast ++) ;
	172
	173	//Now TPC tracks
4b598bb1	174	if(fESD){
	175	//Do it ones, only first time
	176	if(fModule==1){
	177	Int_t nTracks = fESD->GetNumberOfTracks();
c978a69b	178
9a2cdbdf	179	Int_t nPHOSmod = fGeom->GetNModules() ;
c978a69b	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 ;
4b598bb1	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] ;
9a2cdbdf	191	Double_t rEMC = fGeom->GetIPtoCrystalSurface() ; //Use here ideal geometry
4b598bb1	192	for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
4b598bb1	193	track = fESD->GetTrack(iTrack);
c978a69b	194	for(Int_t iTestMod=1; iTestMod<=nPHOSmod; iTestMod++){
9a2cdbdf	195	Double_t modAngle=270.+fGeom->GetPHOSAngle(iTestMod) ;
c978a69b	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 ;
9a2cdbdf	206	fGeom->ImpactOnEmc(vtx, inPHOS.Theta(), inPHOS.Phi(), modNum, z, x) ;
c978a69b	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	}
4b598bb1	215	}
4b598bb1	216	}
861b23c3	217	}
861b23c3	218	}
4b598bb1	219
861b23c3	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
c978a69b	230	TVector3 emcGlobal; // Global position of the EMC recpoint
9a2cdbdf	231	fGeom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);
4b598bb1	232
c978a69b	233	//Calculate actual distance to the PHOS surface
9a2cdbdf	234	Double_t modAngle=270.+fGeom->GetPHOSAngle(emcClu->GetPHOSMod()) ;
c978a69b	235	modAngle*=TMath::Pi()/180. ;
	236	Double_t rEMC = emcGlobal.X()TMath::Cos(modAngle)+emcGlobal.Y()TMath::Sin(modAngle) ;
	237	track->Rotate(modAngle);
861b23c3	238	Double_t xyz[3] ;
c978a69b	239	if(!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)){
	240	dx=999. ;
	241	dz=999. ;
	242	return ; //track coord on the cylinder of PHOS radius
861b23c3	243	}
c978a69b	244	if((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0){
861b23c3	245	dx=999. ;
861b23c3	246	dz=999. ;
c978a69b	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] ;
861b23c3	253
	254	return ;
	255	}
	256	//____________________________________________________________________________
	257	void AliPHOSTrackSegmentMakerv2::Init()
	258	{
	259	// Make all memory allocations that are not possible in default constructor
	260
861b23c3	261	fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
9a2cdbdf	262	fTrackSegments = new TClonesArray("AliPHOSTrackSegment",100);
861b23c3	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 ;
861b23c3	274	}
	275
	276
	277	//____________________________________________________________________________
c978a69b	278	void AliPHOSTrackSegmentMakerv2::MakeLinks()
861b23c3	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
861b23c3	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++ ) {
9a2cdbdf	292	emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP)) ;
c978a69b	293	TVector3 vecEmc ;
c978a69b	294	emcclu->GetLocalPosition(vecEmc) ;
861b23c3	295	Int_t mod=emcclu->GetPHOSMod() ;
c978a69b	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(rawdxrawdx+rawdzrawdz)<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() ;
861b23c3	306	new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(dx, dz, iEmcRP, itrack, -1) ;
	307	}
	308	}
	309	}
861b23c3	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
861b23c3	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
c978a69b	331	Bool_t * tpcExist = 0;
c978a69b	332	Int_t nTracks = fESD->GetNumberOfTracks();
861b23c3	333	if(nTracks>0)
c978a69b	334	tpcExist = new Bool_t[nTracks] ;
861b23c3	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) ;
9a2cdbdf	354	new ((* fTrackSegments)[fNTrackSegments])
9a2cdbdf	355	AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) ,
861b23c3	356	nullpointer,
861b23c3	357	linkUp->GetTrack(),dx,dz) ;
9a2cdbdf	358	(dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
861b23c3	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 ){
9a2cdbdf	372	new ((*fTrackSegments)[fNTrackSegments])
9a2cdbdf	373	AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(iEmcRP+fEmcFirst)),
861b23c3	374	nullpointer) ;
9a2cdbdf	375	(dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
861b23c3	376	fNTrackSegments++;
	377	}
	378	}
	379	}
	380	delete [] emcExist ;
	381	delete [] tpcExist ;
	382	}
	383
	384	//____________________________________________________________________________
9a2cdbdf	385	void AliPHOSTrackSegmentMakerv2::Clusters2TrackSegments(Option_t *option)
861b23c3	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
9a2cdbdf	400	//Make some initializations
	401	fNTrackSegments = 0 ;
	402	fEmcFirst = 0 ;
	403	fEmcLast = 0 ;
861b23c3	404
9a2cdbdf	405	fTrackSegments->Clear();
861b23c3	406
9a2cdbdf	407	// if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
861b23c3	408
9a2cdbdf	409	for(fModule = 1; fModule <= fGeom->GetNModules() ; fModule++ ) {
	410	FillOneModule() ;
	411	MakeLinks() ;
	412	MakePairs() ;
	413	}
861b23c3	414
9a2cdbdf	415	if(strstr(option,"deb"))
9a2cdbdf	416	PrintTrackSegments(option);
861b23c3	417
861b23c3	418	if(strstr(option,"tim")){
861b23c3	419	gBenchmark->Stop("PHOSTSMaker");
9a2cdbdf	420	Info("Exec", "took %f seconds for making TS",
	421	gBenchmark->GetCpuTime("PHOSTSMaker"));
	422	}
861b23c3	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 ") ;
861b23c3	441
861b23c3	442	}
861b23c3	443
861b23c3	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
861b23c3	450	Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
9a2cdbdf	451	printf(" Found %d TrackSegments\n", fTrackSegments->GetEntriesFast() );
861b23c3	452
	453	if(strstr(option,"all")) { // printing found TS
	454	printf("TrackSegment # EMC RP# CPV RP#\n") ;
	455	Int_t index;
9a2cdbdf	456	for (index = 0 ; index <fTrackSegments->GetEntriesFast() ; index++) {
9a2cdbdf	457	AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )fTrackSegments->At(index) ;
861b23c3	458	printf(" %d %d %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetTrackIndex() ) ;
	459	}
	460	}
	461	}