[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(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 ) 
    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() ;    

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