[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 "TTree.h"
#include "TBenchmark.h"

// --- Standard library ---

// --- AliRoot header files ---

#include "AliEMCALTrackSegmentMakerv1.h"
#include "AliEMCALTrackSegment.h"
#include "AliEMCALLink.h"
#include "AliEMCALGetter.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)

}

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

}

//____________________________________________________________________________
Float_t  AliEMCALTrackSegmentMakerv1::HowClose(AliEMCALRecPoint * ec, AliEMCALRecPoint * 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;
  TVector3 vecRP;
  ec->GetGlobalPosition(vecEC);
  rp->GetGlobalPosition(vecRP);
  
  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 ) 
    printf("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());
 
  if ( !gime->TrackSegmentMaker() ) {
    gime->PostTrackSegmentMaker(this);
  }
}

//____________________________________________________________________________
void  AliEMCALTrackSegmentMakerv1::InitParameters()
{
  fClose              = 10e-3 ;   
  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() ;    

  AliEMCALRecPoint * pre ;
  AliEMCALRecPoint * eca ;
  AliEMCALRecPoint * hca ;

  Int_t iPRELink  = 0 ;
  Int_t iHCALink  = 0 ;
    
  Int_t iECARP;
  for(iECARP = 0; iECARP < aECARecPoints->GetEntriesFast(); iECARP++ ) {
    eca = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(iECARP)) ;
    Bool_t toofar = kTRUE ;        
    Int_t iPRERP = 0 ;    
    for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) { 
      pre = dynamic_cast<AliEMCALRecPoint *>(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<AliEMCALRecPoint *>(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 ;
  
  AliEMCALRecPoint * 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<AliEMCALRecPoint *>(aECARecPoints->At(linkPRE->GetECA())) , 
			       dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	fNTrackSegments++ ;
	if (gDebug == 2 ) 
	  printf("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<AliEMCALRecPoint *>(aHCARecPoints->At(linkHCA->GetOther())) ) ;
	  if (gDebug == 2 ) 
	    printf("MakePairs: ECAL section with PRE and HCAL sections") ;
	} 	
	if (!found) {
	  new ((* trackSegments)[fNTrackSegments]) 
	    AliEMCALTrackSegment(dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(linkHCA->GetECA())), null,
				 dynamic_cast<AliEMCALRecPoint *>(aHCARecPoints->At(linkHCA->GetOther()))) ; 
	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	fNTrackSegments++ ;
	if (gDebug == 2 ) 
	  printf("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<AliEMCALRecPoint *>(aECARecPoints->At(iECARP)), null, null) ;
	(dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
	fNTrackSegments++;    
	if( gDebug == 2 ) 
	  printf("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() ;
    
    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");
    printf("Exec: took %f seconds for making TS %f seconds per event", 
	 gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
  }
   Unload();
}

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

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

  printf("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()
{
  // 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() ; 


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