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 ---
36 #include "TBenchmark.h"
38 // --- Standard library ---
40 // --- AliRoot header files ---
42 #include "AliEMCALTrackSegmentMakerv1.h"
43 #include "AliEMCALTrackSegment.h"
44 #include "AliEMCALLink.h"
45 #include "AliEMCALGetter.h"
47 ClassImp( AliEMCALTrackSegmentMakerv1)
50 //____________________________________________________________________________
51 AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1() : AliEMCALTrackSegmentMaker()
53 // default ctor (to be used mainly by Streamer)
56 fDefaultInit = kTRUE ;
59 //____________________________________________________________________________
60 AliEMCALTrackSegmentMakerv1::AliEMCALTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
61 :AliEMCALTrackSegmentMaker(alirunFileName, eventFolderName)
67 fDefaultInit = kFALSE ;
71 //____________________________________________________________________________
72 AliEMCALTrackSegmentMakerv1::~AliEMCALTrackSegmentMakerv1()
75 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
79 //____________________________________________________________________________
80 const TString AliEMCALTrackSegmentMakerv1::BranchName() const
86 //____________________________________________________________________________
87 Float_t AliEMCALTrackSegmentMakerv1::HowClose(AliEMCALRecPoint * ec, AliEMCALRecPoint * rp, Bool_t &toofar)const
89 // Calculates the distance between the EMC RecPoint and the PPSD RecPoint
90 // Clusters are sorted in "rows" and "columns" of width 1 cm
93 Float_t delta = 10. ; // large enough to be ineffective ??!
98 ec->GetGlobalPosition(vecEC);
99 rp->GetGlobalPosition(vecRP);
101 Float_t pro = TMath::Abs(1 - (vecEC * vecRP / ( vecEC.Mag() * vecRP.Mag() ))) ;
111 printf("HowClose: ec = %d, rp = %d pro = %f, toofar=%d", ec->GetIndexInList(), rp->GetIndexInList(), pro, toofar ) ;
116 //____________________________________________________________________________
117 void AliEMCALTrackSegmentMakerv1::Init()
119 // Make all memory allocations that are not possible in default constructor
121 AliEMCALGetter* gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
123 if ( !gime->TrackSegmentMaker() ) {
124 gime->PostTrackSegmentMaker(this);
128 //____________________________________________________________________________
129 void AliEMCALTrackSegmentMakerv1::InitParameters()
132 fTrackSegmentsInRun = 0 ;
136 //____________________________________________________________________________
137 void AliEMCALTrackSegmentMakerv1::MakeLinks()const
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
143 /* AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
144 TObjArray * aECARecPoints = gime->ECARecPoints() ;
145 // TObjArray * aPRERecPoints = gime->PRERecPoints() ;
146 //TObjArray * aHCARecPoints = gime->HCARecPoints() ;
148 fPRELinkArray->Clear() ;
149 fHCALinkArray->Clear() ;
151 AliEMCALRecPoint * pre ;
152 AliEMCALRecPoint * eca ;
153 AliEMCALRecPoint * hca ;
159 for(iECARP = 0; iECARP < aECARecPoints->GetEntriesFast(); iECARP++ ) {
160 eca = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(iECARP)) ;
161 Bool_t toofar = kTRUE ;
163 for(iPRERP = 0; iPRERP < aPRERecPoints->GetEntriesFast(); iPRERP++ ) {
164 pre = dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(iPRERP)) ;
165 Float_t prod = HowClose(eca, pre, toofar) ;
169 new ((*fPRELinkArray)[iPRELink++]) AliEMCALLink(prod, iECARP, iPRERP, 0) ;
174 for(iHCARP = 0; iHCARP < aHCARecPoints->GetEntriesFast(); iHCARP++ ) {
175 hca = dynamic_cast<AliEMCALRecPoint *>(aHCARecPoints->At(iHCARP)) ;
176 Float_t prod = HowClose(eca, hca, toofar) ;
180 new ((*fHCALinkArray)[iHCALink++]) AliEMCALLink(prod, iECARP, iHCARP, 1) ;
185 fPRELinkArray->Sort() ; //first links with largest scalar product
186 fHCALinkArray->Sort() ; //first links with largest scalar product
190 //____________________________________________________________________________
191 void AliEMCALTrackSegmentMakerv1::MakePairs()
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".
198 /*AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
199 TObjArray * aECARecPoints = gime->ECARecPoints() ;
200 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
201 TObjArray * aHCARecPoints = gime->HCARecPoints() ;
202 TClonesArray * trackSegments = gime->TrackSegments() ;
204 //Make arrays to mark clusters already chosen
205 Int_t * ecaExist = 0;
206 Int_t nECA = aECARecPoints->GetEntriesFast() ;
208 ecaExist = new Int_t[nECA] ;
211 for(index = 0; index < nECA; index ++)
212 ecaExist[index] = 1 ;
214 Bool_t * preExist = 0;
215 Int_t nPRE = aPRERecPoints->GetEntriesFast() ;
217 preExist = new Bool_t[nPRE] ;
218 for(index = 0; index < nPRE; index ++)
219 preExist[index] = kTRUE ;
221 Bool_t * hcaExist = 0;
222 Int_t nHCA = aHCARecPoints->GetEntriesFast() ;
224 hcaExist = new Bool_t[nHCA] ;
225 for(index = 0; index < nHCA; index ++)
226 hcaExist[index] = kTRUE ;
228 AliEMCALRecPoint * null = 0 ;
229 // Finds the smallest links and makes pairs of PRE and ECAL clusters with largest scalar product
231 TIter nextPRE(fPRELinkArray) ;
232 AliEMCALLink * linkPRE ;
234 while ( (linkPRE = static_cast<AliEMCALLink *>(nextPRE()) ) ){
236 if(ecaExist[linkPRE->GetECA()] != -1){ //without PRE yet
238 if(preExist[linkPRE->GetOther()]){ // PRE still exist
240 new ((* trackSegments)[fNTrackSegments])
241 AliEMCALTrackSegment(dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(linkPRE->GetECA())) ,
242 dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(linkPRE->GetOther())), null) ;
243 (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
246 printf("MakePairs: ECAL section with PRE section") ;
247 ecaExist[linkPRE->GetECA()] = -1 ; //Mark ecal that pre was found
248 //mark PRE recpoint as already used
249 preExist[linkPRE->GetOther()] = kFALSE ;
250 } //if PRE still exist
254 // Finds the smallest links and makes pairs of HCAL and ECAL clusters with largest scalar product
256 TIter nextHCA(fHCALinkArray) ;
257 AliEMCALLink * linkHCA ;
259 while ( (linkHCA = static_cast<AliEMCALLink *>(nextHCA()) ) ){
261 if(ecaExist[linkHCA->GetECA()] != -2){ //without HCAL yet
263 if(hcaExist[linkHCA->GetOther()]){ // HCAL still exist
264 // search among the already existing track segments
266 Bool_t found = kFALSE ;
267 AliEMCALTrackSegment * ts = 0 ;
268 for ( ii = 0 ; ii < fNTrackSegments ; ii++ ) {
269 ts = dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(ii)) ;
270 if ( ts->GetECAIndex() == linkHCA->GetECA() ) {
276 ts->SetHCARecPoint( dynamic_cast<AliEMCALRecPoint *>(aHCARecPoints->At(linkHCA->GetOther())) ) ;
278 printf("MakePairs: ECAL section with PRE and HCAL sections") ;
281 new ((* trackSegments)[fNTrackSegments])
282 AliEMCALTrackSegment(dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(linkHCA->GetECA())), null,
283 dynamic_cast<AliEMCALRecPoint *>(aHCARecPoints->At(linkHCA->GetOther()))) ;
284 (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
287 printf("MakePairs: ECAL section with HCAL section") ;
289 ecaExist[linkHCA->GetECA()] = -2 ; //Mark ecal that hcal was found
290 //mark HCAL recpoint as already used
291 hcaExist[linkHCA->GetOther()] = kFALSE ;
292 } //if HCAL still exist
297 //look through ECAL recPoints left without PRE/HCAL
298 if(ecaExist){ //if there is ecal rec point
300 for(iECARP = 0; iECARP < nECA ; iECARP++ ){
301 if(ecaExist[iECARP] > 0 ){
302 new ((*trackSegments)[fNTrackSegments])
303 AliEMCALTrackSegment(dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(iECARP)), null, null) ;
304 (dynamic_cast<AliEMCALTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
307 printf("MakePairs: ECAL section alone") ;
317 //____________________________________________________________________________
318 void AliEMCALTrackSegmentMakerv1::Exec(Option_t * option)
323 if(strstr(option,"tim"))
324 gBenchmark->Start("EMCALTSMaker");
326 if(strstr(option,"print")) {
331 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
333 Int_t nevents = gime->MaxEvent() ;
336 for(ievent = 0; ievent < nevents; ievent++){
337 gime->Event(ievent,"R") ;
338 //Make some initializations
339 fNTrackSegments = 0 ;
341 gime->TrackSegments()->Clear() ;
346 WriteTrackSegments() ;
348 if(strstr(option,"deb"))
349 PrintTrackSegments(option) ;
351 //increment the total number of track segments per run
352 fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
356 if(strstr(option,"tim")){
357 gBenchmark->Stop("EMCALTSMaker");
358 printf("Exec: took %f seconds for making TS %f seconds per event",
359 gBenchmark->GetCpuTime("EMCALTSMaker"), gBenchmark->GetCpuTime("EMCALTSMaker")/nevents) ;
364 //____________________________________________________________________________
365 void AliEMCALTrackSegmentMakerv1::Unload()
367 // Unloads the RecPoints and Tracks
368 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
369 gime->EmcalLoader()->UnloadRecPoints() ;
370 gime->EmcalLoader()->UnloadTracks() ;
373 //____________________________________________________________________________
374 void AliEMCALTrackSegmentMakerv1::Print(Option_t * /*option*/)const
376 // Print TrackSegmentMaker parameters
378 printf("Print: TrackSegmentMakerv1 parameters:") ;
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") ;
385 printf("AliEMCALTrackSegmentMakerv1 not initialized ") ;
388 //____________________________________________________________________________
389 void AliEMCALTrackSegmentMakerv1::WriteTrackSegments()
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.
399 AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
401 TClonesArray * trackSegments = gime->TrackSegments() ;
402 trackSegments->Expand(trackSegments->GetEntriesFast()) ;
404 TTree * treeT = gime->TreeT();
407 Int_t bufferSize = 32000 ;
408 TBranch * tsBranch = treeT->Branch("EMCALTS",&trackSegments,bufferSize);
411 gime->WriteTracks("OVERWRITE");
412 gime->WriteTrackSegmentMaker("OVERWRITE");
416 //____________________________________________________________________________
417 void AliEMCALTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
419 // option deb - prints # of found TrackSegments
420 // option deb all - prints as well indexed of found RecParticles assigned to the TS
422 TClonesArray * trackSegments = AliEMCALGetter::Instance()->TrackSegments() ;
425 printf("PrintTrackSegments: Results from TrackSegmentMaker:") ;
426 printf("nevent: %d\n", gAlice->GetEvNumber()) ;
427 printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
429 if(strstr(option,"all")) { // printing found TS
430 printf("TrackSegment# ECAL RP# \n") ;
432 for (index = 0 ; index < fNTrackSegments ; index++) {
433 AliEMCALTrackSegment * ts = (AliEMCALTrackSegment * )trackSegments->At(index) ;
435 ts->GetIndexInList(), ts->GetECAIndex());