Move the tail cancellation from the digitizer to the clusterizer (Christian)
[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 ---
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
47ClassImp( 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)
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}
78
5502f2ea 79//____________________________________________________________________________
80const TString AliEMCALTrackSegmentMakerv1::BranchName() const
81{
88cb7938 82 return GetName() ;
83
5502f2ea 84}
85
86//____________________________________________________________________________
70a93198 87Float_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//____________________________________________________________________________
117void 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//____________________________________________________________________________
129void AliEMCALTrackSegmentMakerv1::InitParameters()
130{
88cb7938 131 fClose = 10e-3 ;
88cb7938 132 fTrackSegmentsInRun = 0 ;
5502f2ea 133}
134
135
136//____________________________________________________________________________
137void 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() ;
146 //TObjArray * aHCARecPoints = gime->HCARecPoints() ;
5502f2ea 147
148 fPRELinkArray->Clear() ;
88cb7938 149 fHCALinkArray->Clear() ;
5502f2ea 150
70a93198 151 AliEMCALRecPoint * pre ;
152 AliEMCALRecPoint * eca ;
153 AliEMCALRecPoint * hca ;
5502f2ea 154
155 Int_t iPRELink = 0 ;
88cb7938 156 Int_t iHCALink = 0 ;
5502f2ea 157
88cb7938 158 Int_t iECARP;
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 ;
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//____________________________________________________________________________
191void 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() ;
199 TObjArray * aECARecPoints = gime->ECARecPoints() ;
5502f2ea 200 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
88cb7938 201 TObjArray * aHCARecPoints = gime->HCARecPoints() ;
202 TClonesArray * trackSegments = gime->TrackSegments() ;
5502f2ea 203
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
210 Int_t index;
88cb7938 211 for(index = 0; index < nECA; index ++)
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())) ,
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) ;
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,
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//____________________________________________________________________________
318void 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
339 fNTrackSegments = 0 ;
5502f2ea 340
88cb7938 341 gime->TrackSegments()->Clear() ;
5502f2ea 342
88cb7938 343 MakeLinks() ;
5502f2ea 344 MakePairs() ;
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) ;
360 }
88cb7938 361 Unload();
362}
363
364//____________________________________________________________________________
365void 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 374void 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 389void 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
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//____________________________________________________________________________
417void 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
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",
435 ts->GetIndexInList(), ts->GetECAIndex());
5502f2ea 436 }
437 }
438}