[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTrackSegmentMakerv1.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$ */
//_________________________________________________________________________
// Implementation version 1 of algorithm class to construct EMCAL track segments
// Track segment for EMCAL is list of 
//        ECAL RecPoint + (possibly) PRE RecPoint + (possibly) HCAL RecPoint
// To find TrackSegments we do the following: 
//  for each ECAL RecPoint we look for PRE and HC RecPoints with same direction  within fSame. 
//  If there is such a PRE or ECAL RecPoint, 
//   we make a "Link": indexes of ECAL and PRE, HCAL  RecPoints and their scalar product. 
//  Then we sort "Links", starting from the 
//   least "Link" pointing to the unassigned RecPoints assigning them to a new TrackSegment. 
//  If there is no PRE, HCAL RecPoint we make a TrackSegment 
//   consisting from ECAL alone. There is no TrackSegments without ECAL RecPoint.
//// In principle this class should be called from AliEMCALReconstructioner, but 
// one can use it as well in standalone mode.
//                 
//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) 
//

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

// --- Standard library ---

// --- AliRoot header files ---

#include "AliEMCALTrackSegmentMakerv1.h"
#include "AliEMCALClusterizerv1.h"
#include "AliEMCALTrackSegment.h"
#include "AliEMCALLink.h"
#include "AliEMCALGetter.h"
#include "AliEMCAL.h"

ClassImp( AliEMCALTrackSegmentMakerv1) 


//____________________________________________________________________________
  AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1() : AliEMCALTrackSegmentMaker()
{
  // default ctor (to be used mainly by Streamer)

  InitParameters() ; 
  fDefaultInit = kTRUE ; 
}

//____________________________________________________________________________
 AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
   :AliEMCALTrackSegmentMaker(alirunFileName, eventFolderName)
{
  // ctor

  InitParameters() ; 
  Init() ;
  fDefaultInit = kFALSE ; 

}

//____________________________________________________________________________
 AliEMCALTrackSegmentMakerv1::~AliEMCALTrackSegmentMakerv1()
{ 
  // dtor
  // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
  
  if (!fDefaultInit) {
    delete fPRELinkArray  ;
    delete fHCALinkArray  ;
  }
}

//____________________________________________________________________________
const TString AliEMCALTrackSegmentMakerv1::BranchName() const 
{  
   return GetName() ;

}

//____________________________________________________________________________
Float_t  AliEMCALTrackSegmentMakerv1::HowClose(AliEMCALTowerRecPoint * ec, AliEMCALTowerRecPoint * rp, Bool_t &toofar)const
{
  // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
  // Clusters are sorted in "rows" and "columns" of width 1 cm

  Float_t r = -1. ;
  Float_t delta = 10. ;  // large enough to be ineffective ??! 

 
  TVector3 vecEC = ec->XYZInAlice() ;
  TVector3 vecRP = rp->XYZInAlice() ;
  
  Float_t pro = TMath::Abs(1 - (vecEC * vecRP / ( vecEC.Mag() * vecRP.Mag() ))) ; 

  if(pro <= delta) { 
    r = pro ;
    toofar = kFALSE ;
  } 
  else 
    toofar = kTRUE ;

  if (gDebug == 2 ) 
    Info("HowClose", "ec = %d, rp = %d pro = %f, toofar=%d", ec->GetIndexInList(), rp->GetIndexInList(), pro, toofar ) ; 
 
  return r ;
}

//____________________________________________________________________________
void  AliEMCALTrackSegmentMakerv1::Init()
{
  // Make all memory allocations that are not possible in default constructor
  
  AliEMCALGetter* gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());

  fPRELinkArray = new TClonesArray("AliEMCALLink", 1000); 
  fHCALinkArray = new TClonesArray("AliEMCALLink", 1000); 
  if ( !gime->TrackSegmentMaker() ) {
    gime->PostTrackSegmentMaker(this);
  }
}

//____________________________________________________________________________
void  AliEMCALTrackSegmentMakerv1::InitParameters()
{
  fClose              = 10e-3 ;   
  fPRELinkArray       = 0 ;
  fHCALinkArray       = 0 ;
  fTrackSegmentsInRun = 0 ; 
}


//____________________________________________________________________________
void  AliEMCALTrackSegmentMakerv1::MakeLinks()const
{ 
  // Finds distances (links) between all PRE, EC and HC clusters, 
  // which are not further apart from each other than fDangle 
  // and sort them in accordance with this distance
  
  AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
  TObjArray * aECARecPoints  = gime->ECARecPoints() ; 
  TObjArray * aPRERecPoints  = gime->PRERecPoints() ; 
  TObjArray * aHCARecPoints  = gime->HCARecPoints() ; 

  fPRELinkArray->Clear() ;    
  fHCALinkArray->Clear() ;    

  AliEMCALTowerRecPoint * pre ;
  AliEMCALTowerRecPoint * eca ;
  AliEMCALTowerRecPoint * hca ;

  Int_t iPRELink  = 0 ;
  Int_t iHCALink  = 0 ;
    
  Int_t iECARP;
  for(iECARP = 0; iECARP < aECARecPoints->GetEntriesFast(); iECARP++ ) {
    eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)) ;
    Bool_t toofar = kTRUE ;        
    Int_t iPRERP = 0 ;    
    for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) { 
      pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(iPRERP)) ;
      Float_t prod = HowClose(eca, pre, toofar) ;    
      if(toofar)
	break ;	 
      if(prod < fClose) { 
	new ((*fPRELinkArray)[iPRELink++])  AliEMCALLink(prod, iECARP, iPRERP, 0) ;
      }      
    }
    toofar = kTRUE ; 
    Int_t iHCARP = 0 ;    
    for(iHCARP = 0; iHCARP < aHCARecPoints->GetEntriesFast(); iHCARP++ ) { 
      hca = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(iHCARP)) ;
      Float_t prod = HowClose(eca, hca, toofar) ;    
      if(toofar)
	break ;	 
      if(prod < fClose) { 
	new ((*fHCALinkArray)[iHCALink++])  AliEMCALLink(prod, iECARP, iHCARP, 1) ;
      }      
    }
  } 
  
  fPRELinkArray->Sort() ;  //first links with largest scalar product
  fHCALinkArray->Sort() ;  //first links with largest scalar product
}

//____________________________________________________________________________
void  AliEMCALTrackSegmentMakerv1::MakePairs()
{ 
  // Using the previously made list of "links", we found the best link - i.e. 
  // link with the largest scalar product (closest to one) to still 
  // unassigned RecParticles. We assign these RecPoints to TrackSegment and 
  // remove them from the list of "unassigned". 
  
  AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 
  TObjArray * aECARecPoints = gime->ECARecPoints() ; 
  TObjArray * aPRERecPoints = gime->PRERecPoints() ; 
  TObjArray * aHCARecPoints = gime->HCARecPoints() ; 
  TClonesArray * trackSegments = gime->TrackSegments() ;   
    
  //Make arrays to mark clusters already chosen
  Int_t * ecaExist = 0;
  Int_t nECA = aECARecPoints->GetEntriesFast() ;  
  if (nECA) 
    ecaExist = new Int_t[nECA] ;
  
  Int_t index;
  for(index = 0; index < nECA; index ++)
    ecaExist[index] = 1 ;
  
  Bool_t * preExist = 0;
  Int_t nPRE = aPRERecPoints->GetEntriesFast() ;  
  if(nPRE)
    preExist = new Bool_t[nPRE] ;
  for(index = 0; index < nPRE; index ++)
    preExist[index] = kTRUE ;
  
  Bool_t * hcaExist = 0;
  Int_t nHCA = aHCARecPoints->GetEntriesFast() ;  
  if(nHCA)
    hcaExist = new Bool_t[nHCA] ;
  for(index = 0; index < nHCA; index ++)
    hcaExist[index] = kTRUE ;
  
  AliEMCALTowerRecPoint * null = 0 ; 
 // Finds the smallest links and makes pairs of PRE and ECAL clusters with largest scalar product 
 
  TIter nextPRE(fPRELinkArray) ;
  AliEMCALLink * linkPRE ;
  
  while ( (linkPRE =  static_cast<AliEMCALLink *>(nextPRE()) ) ){  

    if(ecaExist[linkPRE->GetECA()] != -1){ //without PRE yet 

      if(preExist[linkPRE->GetOther()]){ // PRE still exist
	
	new ((* trackSegments)[fNTrackSegments]) 
	  AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkPRE->GetECA())) , 
			       dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	fNTrackSegments++ ;
	if (gDebug == 2 ) 
	  Info("MakePairs", "ECAL section with PRE section") ; 	
	ecaExist[linkPRE->GetECA()] = -1 ; //Mark ecal  that pre was found 
	//mark PRE recpoint as already used 
	preExist[linkPRE->GetOther()] = kFALSE ;
      } //if PRE still exist
    } 
  }

  // Finds the smallest links and makes pairs of HCAL and ECAL clusters with largest scalar product 
 
  TIter nextHCA(fHCALinkArray) ;
  AliEMCALLink * linkHCA ;
  
  while ( (linkHCA =  static_cast<AliEMCALLink *>(nextHCA()) ) ){  

    if(ecaExist[linkHCA->GetECA()] != -2){ //without HCAL yet 

      if(hcaExist[linkHCA->GetOther()]){ // HCAL still exist
	// search among the already existing track segments
	Int_t ii ; 
	Bool_t found = kFALSE ; 
	AliEMCALTrackSegment * ts = 0 ; 
	for ( ii = 0 ; ii < fNTrackSegments ; ii++ ) {
	  ts = dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(ii)) ;
	  if ( ts->GetECAIndex() == linkHCA->GetECA() ) {
	    found = kTRUE ; 
	    break ; 
	  }
	}
	if (found){
	  ts->SetHCARecPoint( dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther())) ) ;
	  if (gDebug == 2 ) 
	    Info("MakePairs", "ECAL section with PRE and HCAL sections") ;
	} 	
	if (!found) {
	  new ((* trackSegments)[fNTrackSegments]) 
	    AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkHCA->GetECA())), null,
				 dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther()))) ; 
	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	fNTrackSegments++ ;
	if (gDebug == 2 ) 
	  Info("MakePairs", "ECAL section with HCAL section") ; 	
	}
	ecaExist[linkHCA->GetECA()] = -2 ; //Mark ecal  that hcal was found 
	//mark HCAL recpoint as already used 
	hcaExist[linkHCA->GetOther()] = kFALSE ;
      } //if HCAL still exist
    } 
  }
  

  //look through ECAL recPoints left without PRE/HCAL
  if(ecaExist){ //if there is ecal rec point
    Int_t iECARP ;
    for(iECARP = 0; iECARP < nECA  ; iECARP++ ){
      if(ecaExist[iECARP] > 0 ){
	new ((*trackSegments)[fNTrackSegments])  
	  AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)), null, null) ;
	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	fNTrackSegments++;    
	if( gDebug == 2 ) 
	  Info("MakePairs", "ECAL section alone") ; 
     } 
    }
  }
  delete [] ecaExist ; 
  delete [] preExist ; 
  delete [] hcaExist ; 
}

//____________________________________________________________________________
void  AliEMCALTrackSegmentMakerv1::Exec(Option_t * option)
{
  // STEERing method


  if(strstr(option,"tim"))
    gBenchmark->Start("EMCALTSMaker");
 
  if(strstr(option,"print")) {
    Print("") ; 
    return ; 
  }

  AliEMCALGetter * gime = AliEMCALGetter::Instance() ; 

  Int_t nevents = gime->MaxEvent() ;   
  Int_t ievent ;

  for(ievent = 0; ievent < nevents; ievent++){
    gime->Event(ievent,"R") ;
    //Make some initializations 
    fNTrackSegments = 0 ;
    
    gime->TrackSegments()->Clear() ; 
    
    MakeLinks() ;
    MakePairs() ;

    WriteTrackSegments(ievent) ;
    
    if(strstr(option,"deb"))
      PrintTrackSegments(option) ;
    
    //increment the total number of track segments per run 
    fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ; 
    
  }
  
  if(strstr(option,"tim")){
    gBenchmark->Stop("EMCALTSMaker");
    Info("Exec", "took %f seconds for making TS %f seconds per event", 
	 gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
  }
   Unload();
}

//____________________________________________________________________________
void AliEMCALTrackSegmentMakerv1::Unload() 
{
  AliEMCALGetter * gime = AliEMCALGetter::Instance() ;  
  gime->EmcalLoader()->UnloadRecPoints() ;
  gime->EmcalLoader()->UnloadTracks() ;
}

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

  Info("Print", "TrackSegmentMakerv1 parameters:") ; 
  if( strcmp(GetName(), "") != 0 ) { 
    printf("Making Track segments with parameters:\n") ; 
    printf("    Allowed spred on the scalar product of two recpoints with same direction: %f\n", fClose) ;
    printf("============================================\n") ;
  }
  else
    printf("AliEMCALTrackSegmentMakerv1 not initialized ") ;
}

//____________________________________________________________________________
void AliEMCALTrackSegmentMakerv1::WriteTrackSegments(Int_t event)
{
  // Writes found TrackSegments to TreeR. Creates branches 
  // "EMCALTS" and "AliEMCALTrackSegmentMaker" with the same title.
  // In the former branch found TrackSegments are stored, while 
  // in the latter all parameters, with which TS were made. 
  // ROOT does not allow overwriting existing branches, therefore
  // first we check, if branches with the same title already exist.
  // If yes - exits without writing.
  
  AliEMCALGetter *gime = AliEMCALGetter::Instance() ; 

  TClonesArray * trackSegments = gime->TrackSegments() ; 
  trackSegments->Expand(trackSegments->GetEntriesFast()) ;

  TTree * treeT = gime->TreeT();
  
  //First TS
  Int_t bufferSize = 32000 ;    
  TBranch * tsBranch = treeT->Branch("EMCALTS",&trackSegments,bufferSize);
  tsBranch->Fill() ; 

  gime->WriteTracks("OVERWRITE");
  gime->WriteTrackSegmentMaker("OVERWRITE");
}


//____________________________________________________________________________
void AliEMCALTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
{
  // option deb - prints # of found TrackSegments
  // option deb all - prints as well indexed of found RecParticles assigned to the TS
  
  TClonesArray * trackSegments = AliEMCALGetter::Instance()->TrackSegments() ; 


  Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
  printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
  printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 

  if(strstr(option,"all")) {  // printing found TS
    printf("TrackSegment#  ECAL RP#  PRE RP#   HCAL RP#  \n") ; 
    Int_t index;
    for (index = 0 ; index < fNTrackSegments ; index++) {
      AliEMCALTrackSegment * ts = (AliEMCALTrackSegment * )trackSegments->At(index) ; 
      printf("   %d           %d        %d         %d \n", 
	     ts->GetIndexInList(), ts->GetECAIndex(), ts->GetPREIndex(), ts->GetHCAIndex() ); 
    }	
  }
}
Commit	Line	Data
5502f2ea	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	// Implementation version 1 of algorithm class to construct EMCAL track segments
	18	// Track segment for EMCAL is list of
	19	// ECAL RecPoint + (possibly) PRE RecPoint + (possibly) HCAL RecPoint
	20	// To find TrackSegments we do the following:
	21	// for each ECAL RecPoint we look for PRE and HC RecPoints with same direction within fSame.
	22	// If there is such a PRE or ECAL RecPoint,
	23	// we make a "Link": indexes of ECAL and PRE, HCAL RecPoints and their scalar product.
	24	// Then we sort "Links", starting from the
	25	// least "Link" pointing to the unassigned RecPoints assigning them to a new TrackSegment.
	26	// If there is no PRE, HCAL RecPoint we make a TrackSegment
	27	// consisting from ECAL alone. There is no TrackSegments without ECAL RecPoint.
	28	//// In principle this class should be called from AliEMCALReconstructioner, but
	29	// one can use it as well in standalone mode.
	30	//
	31	//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
	32	//
	33
	34	// --- ROOT system ---
	35	#include "TROOT.h"
	36	#include "TFile.h"
	37	#include "TFolder.h"
	38	#include "TTree.h"
	39	#include "TSystem.h"
	40	#include "TBenchmark.h"
	41
	42	// --- Standard library ---
	43
	44	// --- AliRoot header files ---
	45
	46	#include "AliEMCALTrackSegmentMakerv1.h"
	47	#include "AliEMCALClusterizerv1.h"
	48	#include "AliEMCALTrackSegment.h"
	49	#include "AliEMCALLink.h"
	50	#include "AliEMCALGetter.h"
	51	#include "AliEMCAL.h"
5502f2ea	52
	53	ClassImp( AliEMCALTrackSegmentMakerv1)
	54
	55
	56	//____________________________________________________________________________
	57	AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1() : AliEMCALTrackSegmentMaker()
	58	{
	59	// default ctor (to be used mainly by Streamer)
	60
	61	InitParameters() ;
5502f2ea	62	fDefaultInit = kTRUE ;
	63	}
	64
	65	//____________________________________________________________________________
88cb7938	66	AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
88cb7938	67	:AliEMCALTrackSegmentMaker(alirunFileName, eventFolderName)
5502f2ea	68	{
	69	// ctor
	70
	71	InitParameters() ;
5502f2ea	72	Init() ;
5502f2ea	73	fDefaultInit = kFALSE ;
	74
	75	}
	76
	77	//____________________________________________________________________________
	78	AliEMCALTrackSegmentMakerv1::~AliEMCALTrackSegmentMakerv1()
	79	{
	80	// dtor
	81	// fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
	82
	83	if (!fDefaultInit) {
	84	delete fPRELinkArray ;
88cb7938	85	delete fHCALinkArray ;
5502f2ea	86	}
	87	}
	88
5502f2ea	89	//____________________________________________________________________________
	90	const TString AliEMCALTrackSegmentMakerv1::BranchName() const
	91	{
88cb7938	92	return GetName() ;
88cb7938	93
5502f2ea	94	}
	95
	96	//____________________________________________________________________________
	97	Float_t AliEMCALTrackSegmentMakerv1::HowClose(AliEMCALTowerRecPoint * ec, AliEMCALTowerRecPoint * rp, Bool_t &toofar)const
	98	{
	99	// Calculates the distance between the EMC RecPoint and the PPSD RecPoint
	100	// Clusters are sorted in "rows" and "columns" of width 1 cm
	101
	102	Float_t r = -1. ;
	103	Float_t delta = 10. ; // large enough to be ineffective ??!
	104
	105
	106	TVector3 vecEC = ec->XYZInAlice() ;
	107	TVector3 vecRP = rp->XYZInAlice() ;
	108
	109	Float_t pro = TMath::Abs(1 - (vecEC * vecRP / ( vecEC.Mag() * vecRP.Mag() ))) ;
	110
	111	if(pro <= delta) {
	112	r = pro ;
	113	toofar = kFALSE ;
	114	}
	115	else
	116	toofar = kTRUE ;
	117
	118	if (gDebug == 2 )
	119	Info("HowClose", "ec = %d, rp = %d pro = %f, toofar=%d", ec->GetIndexInList(), rp->GetIndexInList(), pro, toofar ) ;
	120
	121	return r ;
	122	}
	123
	124	//____________________________________________________________________________
	125	void AliEMCALTrackSegmentMakerv1::Init()
	126	{
	127	// Make all memory allocations that are not possible in default constructor
	128
88cb7938	129	AliEMCALGetter* gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
5502f2ea	130
88cb7938	131	fPRELinkArray = new TClonesArray("AliEMCALLink", 1000);
	132	fHCALinkArray = new TClonesArray("AliEMCALLink", 1000);
	133	if ( !gime->TrackSegmentMaker() ) {
	134	gime->PostTrackSegmentMaker(this);
	135	}
5502f2ea	136	}
	137
	138	//____________________________________________________________________________
	139	void AliEMCALTrackSegmentMakerv1::InitParameters()
	140	{
88cb7938	141	fClose = 10e-3 ;
	142	fPRELinkArray = 0 ;
	143	fHCALinkArray = 0 ;
	144	fTrackSegmentsInRun = 0 ;
5502f2ea	145	}
	146
	147
	148	//____________________________________________________________________________
	149	void AliEMCALTrackSegmentMakerv1::MakeLinks()const
	150	{
	151	// Finds distances (links) between all PRE, EC and HC clusters,
	152	// which are not further apart from each other than fDangle
	153	// and sort them in accordance with this distance
	154
88cb7938	155	AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
	156	TObjArray * aECARecPoints = gime->ECARecPoints() ;
	157	TObjArray * aPRERecPoints = gime->PRERecPoints() ;
	158	TObjArray * aHCARecPoints = gime->HCARecPoints() ;
5502f2ea	159
5502f2ea	160	fPRELinkArray->Clear() ;
88cb7938	161	fHCALinkArray->Clear() ;
5502f2ea	162
5502f2ea	163	AliEMCALTowerRecPoint * pre ;
88cb7938	164	AliEMCALTowerRecPoint * eca ;
88cb7938	165	AliEMCALTowerRecPoint * hca ;
5502f2ea	166
5502f2ea	167	Int_t iPRELink = 0 ;
88cb7938	168	Int_t iHCALink = 0 ;
5502f2ea	169
88cb7938	170	Int_t iECARP;
	171	for(iECARP = 0; iECARP < aECARecPoints->GetEntriesFast(); iECARP++ ) {
	172	eca = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)) ;
5502f2ea	173	Bool_t toofar = kTRUE ;
	174	Int_t iPRERP = 0 ;
	175	for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) {
	176	pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(iPRERP)) ;
88cb7938	177	Float_t prod = HowClose(eca, pre, toofar) ;
5502f2ea	178	if(toofar)
	179	break ;
	180	if(prod < fClose) {
88cb7938	181	new ((*fPRELinkArray)[iPRELink++]) AliEMCALLink(prod, iECARP, iPRERP, 0) ;
5502f2ea	182	}
	183	}
	184	toofar = kTRUE ;
88cb7938	185	Int_t iHCARP = 0 ;
	186	for(iHCARP = 0; iHCARP < aHCARecPoints->GetEntriesFast(); iHCARP++ ) {
	187	hca = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(iHCARP)) ;
	188	Float_t prod = HowClose(eca, hca, toofar) ;
5502f2ea	189	if(toofar)
	190	break ;
	191	if(prod < fClose) {
88cb7938	192	new ((*fHCALinkArray)[iHCALink++]) AliEMCALLink(prod, iECARP, iHCARP, 1) ;
5502f2ea	193	}
	194	}
	195	}
	196
	197	fPRELinkArray->Sort() ; //first links with largest scalar product
88cb7938	198	fHCALinkArray->Sort() ; //first links with largest scalar product
5502f2ea	199	}
	200
	201	//____________________________________________________________________________
	202	void AliEMCALTrackSegmentMakerv1::MakePairs()
	203	{
	204	// Using the previously made list of "links", we found the best link - i.e.
	205	// link with the largest scalar product (closest to one) to still
	206	// unassigned RecParticles. We assign these RecPoints to TrackSegment and
	207	// remove them from the list of "unassigned".
	208
88cb7938	209	AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
88cb7938	210	TObjArray * aECARecPoints = gime->ECARecPoints() ;
5502f2ea	211	TObjArray * aPRERecPoints = gime->PRERecPoints() ;
88cb7938	212	TObjArray * aHCARecPoints = gime->HCARecPoints() ;
88cb7938	213	TClonesArray * trackSegments = gime->TrackSegments() ;
5502f2ea	214
5502f2ea	215	//Make arrays to mark clusters already chosen
88cb7938	216	Int_t * ecaExist = 0;
	217	Int_t nECA = aECARecPoints->GetEntriesFast() ;
	218	if (nECA)
	219	ecaExist = new Int_t[nECA] ;
5502f2ea	220
5502f2ea	221	Int_t index;
88cb7938	222	for(index = 0; index < nECA; index ++)
88cb7938	223	ecaExist[index] = 1 ;
5502f2ea	224
	225	Bool_t * preExist = 0;
	226	Int_t nPRE = aPRERecPoints->GetEntriesFast() ;
	227	if(nPRE)
	228	preExist = new Bool_t[nPRE] ;
	229	for(index = 0; index < nPRE; index ++)
	230	preExist[index] = kTRUE ;
	231
88cb7938	232	Bool_t * hcaExist = 0;
	233	Int_t nHCA = aHCARecPoints->GetEntriesFast() ;
	234	if(nHCA)
	235	hcaExist = new Bool_t[nHCA] ;
	236	for(index = 0; index < nHCA; index ++)
	237	hcaExist[index] = kTRUE ;
5502f2ea	238
	239	AliEMCALTowerRecPoint * null = 0 ;
	240	// Finds the smallest links and makes pairs of PRE and ECAL clusters with largest scalar product
	241
	242	TIter nextPRE(fPRELinkArray) ;
	243	AliEMCALLink * linkPRE ;
	244
	245	while ( (linkPRE = static_cast<AliEMCALLink *>(nextPRE()) ) ){
	246
88cb7938	247	if(ecaExist[linkPRE->GetECA()] != -1){ //without PRE yet
5502f2ea	248
	249	if(preExist[linkPRE->GetOther()]){ // PRE still exist
	250
	251	new ((* trackSegments)[fNTrackSegments])
88cb7938	252	AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkPRE->GetECA())) ,
5502f2ea	253	dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
	254	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	255	fNTrackSegments++ ;
	256	if (gDebug == 2 )
	257	Info("MakePairs", "ECAL section with PRE section") ;
88cb7938	258	ecaExist[linkPRE->GetECA()] = -1 ; //Mark ecal that pre was found
5502f2ea	259	//mark PRE recpoint as already used
	260	preExist[linkPRE->GetOther()] = kFALSE ;
	261	} //if PRE still exist
	262	}
	263	}
	264
	265	// Finds the smallest links and makes pairs of HCAL and ECAL clusters with largest scalar product
	266
88cb7938	267	TIter nextHCA(fHCALinkArray) ;
88cb7938	268	AliEMCALLink * linkHCA ;
5502f2ea	269
88cb7938	270	while ( (linkHCA = static_cast<AliEMCALLink *>(nextHCA()) ) ){
5502f2ea	271
88cb7938	272	if(ecaExist[linkHCA->GetECA()] != -2){ //without HCAL yet
5502f2ea	273
88cb7938	274	if(hcaExist[linkHCA->GetOther()]){ // HCAL still exist
5502f2ea	275	// search among the already existing track segments
	276	Int_t ii ;
	277	Bool_t found = kFALSE ;
	278	AliEMCALTrackSegment * ts = 0 ;
	279	for ( ii = 0 ; ii < fNTrackSegments ; ii++ ) {
	280	ts = dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(ii)) ;
88cb7938	281	if ( ts->GetECAIndex() == linkHCA->GetECA() ) {
5502f2ea	282	found = kTRUE ;
	283	break ;
	284	}
	285	}
	286	if (found){
88cb7938	287	ts->SetHCARecPoint( dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther())) ) ;
5502f2ea	288	if (gDebug == 2 )
	289	Info("MakePairs", "ECAL section with PRE and HCAL sections") ;
	290	}
	291	if (!found) {
	292	new ((* trackSegments)[fNTrackSegments])
88cb7938	293	AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(linkHCA->GetECA())), null,
88cb7938	294	dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(linkHCA->GetOther()))) ;
5502f2ea	295	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	296	fNTrackSegments++ ;
	297	if (gDebug == 2 )
	298	Info("MakePairs", "ECAL section with HCAL section") ;
	299	}
88cb7938	300	ecaExist[linkHCA->GetECA()] = -2 ; //Mark ecal that hcal was found
5502f2ea	301	//mark HCAL recpoint as already used
88cb7938	302	hcaExist[linkHCA->GetOther()] = kFALSE ;
5502f2ea	303	} //if HCAL still exist
	304	}
	305	}
	306
	307
	308	//look through ECAL recPoints left without PRE/HCAL
88cb7938	309	if(ecaExist){ //if there is ecal rec point
	310	Int_t iECARP ;
	311	for(iECARP = 0; iECARP < nECA ; iECARP++ ){
	312	if(ecaExist[iECARP] > 0 ){
5502f2ea	313	new ((*trackSegments)[fNTrackSegments])
88cb7938	314	AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(iECARP)), null, null) ;
5502f2ea	315	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	316	fNTrackSegments++;
	317	if( gDebug == 2 )
	318	Info("MakePairs", "ECAL section alone") ;
	319	}
	320	}
	321	}
88cb7938	322	delete [] ecaExist ;
5502f2ea	323	delete [] preExist ;
88cb7938	324	delete [] hcaExist ;
5502f2ea	325	}
	326
	327	//____________________________________________________________________________
	328	void AliEMCALTrackSegmentMakerv1::Exec(Option_t * option)
	329	{
	330	// STEERing method
	331
5502f2ea	332
	333	if(strstr(option,"tim"))
	334	gBenchmark->Start("EMCALTSMaker");
	335
	336	if(strstr(option,"print")) {
	337	Print("") ;
	338	return ;
	339	}
	340
88cb7938	341	AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
5502f2ea	342
88cb7938	343	Int_t nevents = gime->MaxEvent() ;
5502f2ea	344	Int_t ievent ;
	345
	346	for(ievent = 0; ievent < nevents; ievent++){
5502f2ea	347	gime->Event(ievent,"R") ;
5502f2ea	348	//Make some initializations
5502f2ea	349	fNTrackSegments = 0 ;
5502f2ea	350
88cb7938	351	gime->TrackSegments()->Clear() ;
5502f2ea	352
88cb7938	353	MakeLinks() ;
5502f2ea	354	MakePairs() ;
	355
	356	WriteTrackSegments(ievent) ;
	357
	358	if(strstr(option,"deb"))
	359	PrintTrackSegments(option) ;
	360
	361	//increment the total number of track segments per run
88cb7938	362	fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
5502f2ea	363
	364	}
	365
	366	if(strstr(option,"tim")){
	367	gBenchmark->Stop("EMCALTSMaker");
	368	Info("Exec", "took %f seconds for making TS %f seconds per event",
	369	gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
	370	}
88cb7938	371	Unload();
	372	}
	373
	374	//____________________________________________________________________________
	375	void AliEMCALTrackSegmentMakerv1::Unload()
	376	{
	377	AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
	378	gime->EmcalLoader()->UnloadRecPoints() ;
	379	gime->EmcalLoader()->UnloadTracks() ;
5502f2ea	380	}
	381
	382	//____________________________________________________________________________
	383	void AliEMCALTrackSegmentMakerv1::Print(Option_t * option)const
	384	{
	385	// Print TrackSegmentMaker parameters
	386
	387	Info("Print", "TrackSegmentMakerv1 parameters:") ;
	388	if( strcmp(GetName(), "") != 0 ) {
	389	printf("Making Track segments with parameters:\n") ;
	390	printf(" Allowed spred on the scalar product of two recpoints with same direction: %f\n", fClose) ;
	391	printf("============================================\n") ;
	392	}
	393	else
	394	printf("AliEMCALTrackSegmentMakerv1 not initialized ") ;
	395	}
	396
	397	//____________________________________________________________________________
	398	void AliEMCALTrackSegmentMakerv1::WriteTrackSegments(Int_t event)
	399	{
	400	// Writes found TrackSegments to TreeR. Creates branches
	401	// "EMCALTS" and "AliEMCALTrackSegmentMaker" with the same title.
	402	// In the former branch found TrackSegments are stored, while
	403	// in the latter all parameters, with which TS were made.
	404	// ROOT does not allow overwriting existing branches, therefore
	405	// first we check, if branches with the same title already exist.
	406	// If yes - exits without writing.
	407
88cb7938	408	AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
5502f2ea	409
	410	TClonesArray * trackSegments = gime->TrackSegments() ;
	411	trackSegments->Expand(trackSegments->GetEntriesFast()) ;
88cb7938	412
88cb7938	413	TTree * treeT = gime->TreeT();
5502f2ea	414
	415	//First TS
	416	Int_t bufferSize = 32000 ;
88cb7938	417	TBranch * tsBranch = treeT->Branch("EMCALTS",&trackSegments,bufferSize);
	418	tsBranch->Fill() ;
	419
	420	gime->WriteTracks("OVERWRITE");
	421	gime->WriteTrackSegmentMaker("OVERWRITE");
5502f2ea	422	}
	423
	424
	425	//____________________________________________________________________________
	426	void AliEMCALTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
	427	{
	428	// option deb - prints # of found TrackSegments
	429	// option deb all - prints as well indexed of found RecParticles assigned to the TS
5502f2ea	430
88cb7938	431	TClonesArray * trackSegments = AliEMCALGetter::Instance()->TrackSegments() ;
5502f2ea	432
	433
	434	Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
	435	printf("nevent: %d\n", gAlice->GetEvNumber()) ;
	436	printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
	437
	438	if(strstr(option,"all")) { // printing found TS
	439	printf("TrackSegment# ECAL RP# PRE RP# HCAL RP# \n") ;
	440	Int_t index;
	441	for (index = 0 ; index < fNTrackSegments ; index++) {
	442	AliEMCALTrackSegment * ts = (AliEMCALTrackSegment * )trackSegments->At(index) ;
	443	printf(" %d %d %d %d \n",
88cb7938	444	ts->GetIndexInList(), ts->GetECAIndex(), ts->GetPREIndex(), ts->GetHCAIndex() );
5502f2ea	445	}
	446	}
	447	}