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