1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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.
31 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
34 // --- ROOT system ---
40 #include "TBenchmark.h"
42 // --- Standard library ---
44 // --- AliRoot header files ---
46 #include "AliEMCALTrackSegmentMakerv1.h"
47 #include "AliEMCALClusterizerv1.h"
48 #include "AliEMCALTrackSegment.h"
49 #include "AliEMCALLink.h"
50 #include "AliEMCALGetter.h"
54 ClassImp( AliEMCALTrackSegmentMakerv1)
57 //____________________________________________________________________________
58 AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1() : AliEMCALTrackSegmentMaker()
60 // default ctor (to be used mainly by Streamer)
64 fTrackSegmentsInRun = 0 ;
66 fDefaultInit = kTRUE ;
69 //____________________________________________________________________________
70 AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1(const char * headerFile, const char * name, const Bool_t toSplit) : AliEMCALTrackSegmentMaker(headerFile, name, toSplit)
76 fTrackSegmentsInRun = 0 ;
80 fDefaultInit = kFALSE ;
84 //____________________________________________________________________________
85 AliEMCALTrackSegmentMakerv1::~AliEMCALTrackSegmentMakerv1()
88 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
91 delete fPRELinkArray ;
98 //____________________________________________________________________________
99 const TString AliEMCALTrackSegmentMakerv1::BranchName() const
101 TString branchName(GetName() ) ;
102 branchName.Remove(branchName.Index(Version())-1) ;
106 //____________________________________________________________________________
107 Float_t AliEMCALTrackSegmentMakerv1::HowClose(AliEMCALTowerRecPoint * ec, AliEMCALTowerRecPoint * rp, Bool_t &toofar)const
109 // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
110 // Clusters are sorted in "rows" and "columns" of width 1 cm
113 Float_t delta = 10. ; // large enough to be ineffective ??!
116 TVector3 vecEC = ec->XYZInAlice() ;
117 TVector3 vecRP = rp->XYZInAlice() ;
119 Float_t pro = TMath::Abs(1 - (vecEC * vecRP / ( vecEC.Mag() * vecRP.Mag() ))) ;
129 Info("HowClose", "ec = %d, rp = %d pro = %f, toofar=%d", ec->GetIndexInList(), rp->GetIndexInList(), pro, toofar ) ;
134 //____________________________________________________________________________
135 void AliEMCALTrackSegmentMakerv1::Init()
137 // Make all memory allocations that are not possible in default constructor
139 if ( strcmp(GetTitle(), "") == 0 )
140 SetTitle("galice.root") ;
142 TString branchname = GetName() ;
143 branchname.Remove(branchname.Index(Version())-1) ;
144 AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(),branchname.Data(), fToSplit ) ;
146 Error("Init", "Could not obtain the Getter object !") ;
152 //First - extract full path if necessary
153 TString fileName(GetTitle()) ;
154 Ssiz_t islash = fileName.Last('/') ;
155 if(islash<fileName.Length())
156 fileName.Remove(islash+1,fileName.Length()) ;
159 fileName+="EMCAL.RecData." ;
160 if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
161 fileName+=branchname ;
165 fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
167 fSplitFile = TFile::Open(fileName.Data(),"update") ;
170 fPRELinkArray = new TClonesArray("AliEMCALLink", 1000);
171 fHCLinkArray = new TClonesArray("AliEMCALLink", 1000);
174 gime->PostTrackSegmentMaker(this) ;
175 gime->PostTrackSegments(BranchName()) ;
179 //____________________________________________________________________________
180 void AliEMCALTrackSegmentMakerv1::InitParameters()
185 TString tsmName( GetName()) ;
186 if (tsmName.IsNull() )
187 tsmName = "Default" ;
188 tsmName.Append(":") ;
189 tsmName.Append(Version()) ;
194 //____________________________________________________________________________
195 void AliEMCALTrackSegmentMakerv1::MakeLinks()const
197 // Finds distances (links) between all PRE, EC and HC clusters,
198 // which are not further apart from each other than fDangle
199 // and sort them in accordance with this distance
201 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
202 TObjArray * aECRecPoints = gime->ECALRecPoints() ;
203 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
204 TObjArray * aHCRecPoints = gime->HCALRecPoints() ;
206 fPRELinkArray->Clear() ;
207 fHCLinkArray->Clear() ;
209 AliEMCALTowerRecPoint * pre ;
210 AliEMCALTowerRecPoint * ec ;
211 AliEMCALTowerRecPoint * hc ;
217 for(iECRP = 0; iECRP < aECRecPoints->GetEntriesFast(); iECRP++ ) {
218 ec = dynamic_cast<AliEMCALTowerRecPoint *>(aECRecPoints->At(iECRP)) ;
219 Bool_t toofar = kTRUE ;
221 for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) {
222 pre = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(iPRERP)) ;
223 Float_t prod = HowClose(ec, pre, toofar) ;
227 new ((*fPRELinkArray)[iPRELink++]) AliEMCALLink(prod, iECRP, iPRERP, 0) ;
232 for(iHCRP = 0; iHCRP < aHCRecPoints->GetEntriesFast(); iHCRP++ ) {
233 hc = dynamic_cast<AliEMCALTowerRecPoint *>(aHCRecPoints->At(iHCRP)) ;
234 Float_t prod = HowClose(ec, hc, toofar) ;
238 new ((*fHCLinkArray)[iHCLink++]) AliEMCALLink(prod, iECRP, iHCRP, 1) ;
243 fPRELinkArray->Sort() ; //first links with largest scalar product
244 fHCLinkArray->Sort() ; //first links with largest scalar product
247 //____________________________________________________________________________
248 void AliEMCALTrackSegmentMakerv1::MakePairs()
250 // Using the previously made list of "links", we found the best link - i.e.
251 // link with the largest scalar product (closest to one) to still
252 // unassigned RecParticles. We assign these RecPoints to TrackSegment and
253 // remove them from the list of "unassigned".
255 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
256 TObjArray * aECALRecPoints = gime->ECALRecPoints() ;
257 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
258 TObjArray * aHCALRecPoints = gime->HCALRecPoints() ;
259 TClonesArray * trackSegments = gime->TrackSegments(BranchName()) ;
261 //Make arrays to mark clusters already chosen
262 Int_t * ecalExist = 0;
263 Int_t nECAL = aECALRecPoints->GetEntriesFast() ;
265 ecalExist = new Int_t[nECAL] ;
268 for(index = 0; index < nECAL; index ++)
269 ecalExist[index] = 1 ;
271 Bool_t * preExist = 0;
272 Int_t nPRE = aPRERecPoints->GetEntriesFast() ;
274 preExist = new Bool_t[nPRE] ;
275 for(index = 0; index < nPRE; index ++)
276 preExist[index] = kTRUE ;
278 Bool_t * hcalExist = 0;
279 Int_t nHCAL = aHCALRecPoints->GetEntriesFast() ;
281 hcalExist = new Bool_t[nHCAL] ;
282 for(index = 0; index < nHCAL; index ++)
283 hcalExist[index] = kTRUE ;
285 AliEMCALTowerRecPoint * null = 0 ;
286 // Finds the smallest links and makes pairs of PRE and ECAL clusters with largest scalar product
288 TIter nextPRE(fPRELinkArray) ;
289 AliEMCALLink * linkPRE ;
291 while ( (linkPRE = static_cast<AliEMCALLink *>(nextPRE()) ) ){
293 if(ecalExist[linkPRE->GetECAL()] != -1){ //without PRE yet
295 if(preExist[linkPRE->GetOther()]){ // PRE still exist
297 new ((* trackSegments)[fNTrackSegments])
298 AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECALRecPoints->At(linkPRE->GetECAL())) ,
299 dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
300 (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
303 Info("MakePairs", "ECAL section with PRE section") ;
304 ecalExist[linkPRE->GetECAL()] = -1 ; //Mark ecal that pre was found
305 //mark PRE recpoint as already used
306 preExist[linkPRE->GetOther()] = kFALSE ;
307 } //if PRE still exist
311 // Finds the smallest links and makes pairs of HCAL and ECAL clusters with largest scalar product
313 TIter nextHCAL(fHCLinkArray) ;
314 AliEMCALLink * linkHCAL ;
316 while ( (linkHCAL = static_cast<AliEMCALLink *>(nextHCAL()) ) ){
318 if(ecalExist[linkHCAL->GetECAL()] != -2){ //without HCAL yet
320 if(hcalExist[linkHCAL->GetOther()]){ // HCAL still exist
321 // search among the already existing track segments
323 Bool_t found = kFALSE ;
324 AliEMCALTrackSegment * ts = 0 ;
325 for ( ii = 0 ; ii < fNTrackSegments ; ii++ ) {
326 ts = dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(ii)) ;
327 if ( ts->GetECIndex() == linkHCAL->GetECAL() ) {
333 ts->SetHCRecPoint( dynamic_cast<AliEMCALTowerRecPoint *>(aHCALRecPoints->At(linkHCAL->GetOther())) ) ;
335 Info("MakePairs", "ECAL section with PRE and HCAL sections") ;
338 new ((* trackSegments)[fNTrackSegments])
339 AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECALRecPoints->At(linkHCAL->GetECAL())), null,
340 dynamic_cast<AliEMCALTowerRecPoint *>(aHCALRecPoints->At(linkHCAL->GetOther()))) ;
341 (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
344 Info("MakePairs", "ECAL section with HCAL section") ;
346 ecalExist[linkHCAL->GetECAL()] = -2 ; //Mark ecal that hcal was found
347 //mark HCAL recpoint as already used
348 hcalExist[linkHCAL->GetOther()] = kFALSE ;
349 } //if HCAL still exist
354 //look through ECAL recPoints left without PRE/HCAL
355 if(ecalExist){ //if there is ecal rec point
357 for(iECALRP = 0; iECALRP < nECAL ; iECALRP++ ){
358 if(ecalExist[iECALRP] > 0 ){
359 new ((*trackSegments)[fNTrackSegments])
360 AliEMCALTrackSegment(dynamic_cast<AliEMCALTowerRecPoint *>(aECALRecPoints->At(iECALRP)), null, null) ;
361 (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
364 Info("MakePairs", "ECAL section alone") ;
368 delete [] ecalExist ;
370 delete [] hcalExist ;
373 //____________________________________________________________________________
374 void AliEMCALTrackSegmentMakerv1::Exec(Option_t * option)
378 if( strcmp(GetName(), "")== 0 )
381 if(strstr(option,"tim"))
382 gBenchmark->Start("EMCALTSMaker");
384 if(strstr(option,"print")) {
389 AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
390 if(gime->BranchExists("TrackSegments") )
393 Int_t nevents = gime->MaxEvent() ; //(Int_t) gAlice->TreeE()->GetEntries() ;
396 for(ievent = 0; ievent < nevents; ievent++){
398 gime->Event(ievent,"R") ;
400 //Make some initializations
401 fNTrackSegments = 0 ;
402 gime->TrackSegments(BranchName())->Clear() ;
408 WriteTrackSegments(ievent) ;
410 if(strstr(option,"deb"))
411 PrintTrackSegments(option) ;
413 //increment the total number of track segments per run
414 fTrackSegmentsInRun += gime->TrackSegments(BranchName())->GetEntriesFast() ;
418 if(strstr(option,"tim")){
419 gBenchmark->Stop("EMCALTSMaker");
420 Info("Exec", "took %f seconds for making TS %f seconds per event",
421 gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
426 //____________________________________________________________________________
427 void AliEMCALTrackSegmentMakerv1::Print(Option_t * option)const
429 // Print TrackSegmentMaker parameters
431 Info("Print", "TrackSegmentMakerv1 parameters:") ;
432 if( strcmp(GetName(), "") != 0 ) {
433 printf("Making Track segments with parameters:\n") ;
434 printf(" Allowed spred on the scalar product of two recpoints with same direction: %f\n", fClose) ;
435 printf("============================================\n") ;
438 printf("AliEMCALTrackSegmentMakerv1 not initialized ") ;
441 //____________________________________________________________________________
442 void AliEMCALTrackSegmentMakerv1::WriteTrackSegments(Int_t event)
444 // Writes found TrackSegments to TreeR. Creates branches
445 // "EMCALTS" and "AliEMCALTrackSegmentMaker" with the same title.
446 // In the former branch found TrackSegments are stored, while
447 // in the latter all parameters, with which TS were made.
448 // ROOT does not allow overwriting existing branches, therefore
449 // first we check, if branches with the same title already exist.
450 // If yes - exits without writing.
452 AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ;
454 TClonesArray * trackSegments = gime->TrackSegments() ;
455 trackSegments->Expand(trackSegments->GetEntriesFast()) ;
463 sprintf(name,"%s%d", "TreeR",event) ;
464 treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
467 treeR = gAlice->TreeR();
471 gAlice->MakeTree("R", fSplitFile);
472 treeR = gAlice->TreeR() ;
476 Int_t bufferSize = 32000 ;
477 TBranch * tsBranch = treeR->Branch("EMCALTS",&trackSegments,bufferSize);
478 tsBranch->SetTitle(BranchName());
481 Int_t splitlevel = 0 ;
482 AliEMCALTrackSegmentMakerv1 * ts = this ;
483 TBranch * tsMakerBranch = treeR->Branch("AliEMCALTrackSegmentMaker","AliEMCALTrackSegmentMakerv1",
484 &ts,bufferSize,splitlevel);
485 tsMakerBranch->SetTitle(BranchName());
488 tsMakerBranch->Fill() ;
490 treeR->AutoSave() ; //Write(0,kOverwrite) ;
491 if(gAlice->TreeR()!=treeR)
496 //____________________________________________________________________________
497 void AliEMCALTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
499 // option deb - prints # of found TrackSegments
500 // option deb all - prints as well indexed of found RecParticles assigned to the TS
501 TString taskName(GetName()) ;
502 taskName.Remove(taskName.Index(Version())-1) ;
504 TClonesArray * trackSegments = AliEMCALGetter::GetInstance()->TrackSegments(taskName) ;
507 Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
508 printf("nevent: %d\n", gAlice->GetEvNumber()) ;
509 printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
511 if(strstr(option,"all")) { // printing found TS
512 printf("TrackSegment# ECAL RP# PRE RP# HCAL RP# \n") ;
514 for (index = 0 ; index < fNTrackSegments ; index++) {
515 AliEMCALTrackSegment * ts = (AliEMCALTrackSegment * )trackSegments->At(index) ;
516 printf(" %d %d %d %d \n",
517 ts->GetIndexInList(), ts->GetECIndex(), ts->GetPREIndex(), ts->GetHCIndex() );