Memory leak fixed
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv2.cxx
CommitLineData
861b23c3 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/* History of cvs commits:
18 *
19 * $Log$
9ee9f78d 20 * Revision 1.3 2007/04/25 19:39:42 kharlov
21 * Track extracpolation improved
22 *
c978a69b 23 * Revision 1.2 2007/04/01 19:16:52 kharlov
24 * D.P.: Produce EMCTrackSegments using TPC/ITS tracks (no CPV)
25 *
861b23c3 26 *
27 */
28
29//_________________________________________________________________________
30// Implementation version 2 of algorithm class to construct PHOS track segments
31// Track segment for PHOS is list of
32// EMC RecPoint + (possibly) projection of TPC track
33// To find TrackSegments we do the following:
34// for each EMC RecPoints we look at
35// TPC projections radius fRtpc.
36// If there is such a track
37// we make "Link" it is just indexes of EMC and TPC track and distance
38// between them in the PHOS plane.
39// Then we sort "Links" and starting from the
40// least "Link" pointing to the unassined EMC and TPC assing them to
41// new TrackSegment.
42// If there is no TPC track we make TrackSegment
43// consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
44//// In principle this class should be called from AliPHOSReconstructor, but
45// one can use it as well in standalone mode.
46// Use case:
47// root [0] AliPHOSTrackSegmentMakerv2 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
48// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
49// // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
50// // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")
51// root [1] t->ExecuteTask()
52// root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
53// root [4] t->ExecuteTask("deb all time")
54//
55//*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
56//
57
58// --- ROOT system ---
4b598bb1 59#include "TFile.h"
861b23c3 60#include "TTree.h"
61#include "TBenchmark.h"
62
63// --- Standard library ---
64#include "Riostream.h"
65// --- AliRoot header files ---
66#include "AliPHOSGeometry.h"
67#include "AliPHOSTrackSegmentMakerv2.h"
68#include "AliPHOSTrackSegment.h"
69#include "AliPHOSLink.h"
70#include "AliPHOSGetter.h"
71#include "AliESD.h"
72#include "AliESDtrack.h"
73
74ClassImp( AliPHOSTrackSegmentMakerv2)
75
76
77//____________________________________________________________________________
78AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2() :
79 AliPHOSTrackSegmentMaker(),
80 fDefaultInit(kTRUE),
81 fWrite(kFALSE),
82 fNTrackSegments(0),
83 fRtpc(0.f),
c978a69b 84 fVtx(0.,0.,0.),
861b23c3 85 fLinkUpArray(0),
c978a69b 86 fTPCtracks(),
861b23c3 87 fEmcFirst(0),
88 fEmcLast(0),
89 fModule(0),
90 fTrackSegmentsInRun(0)
861b23c3 91{
92 // default ctor (to be used mainly by Streamer)
93 InitParameters() ;
861b23c3 94}
95
96//____________________________________________________________________________
97AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const TString & alirunFileName, const TString & eventFolderName) :
98 AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName),
99 fDefaultInit(kFALSE),
100 fWrite(kFALSE),
101 fNTrackSegments(0),
102 fRtpc(0.f),
c978a69b 103 fVtx(0.,0.,0.),
861b23c3 104 fLinkUpArray(0),
c978a69b 105 fTPCtracks(),
861b23c3 106 fEmcFirst(0),
107 fEmcLast(0),
108 fModule(0),
109 fTrackSegmentsInRun(0)
110{
111 // ctor
112 InitParameters() ;
113 Init() ;
114 fESD = 0;
115}
116
117//____________________________________________________________________________
118AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
119 AliPHOSTrackSegmentMaker(tsm),
120 fDefaultInit(kFALSE),
121 fWrite(kFALSE),
122 fNTrackSegments(0),
123 fRtpc(0.f),
c978a69b 124 fVtx(0.,0.,0.),
861b23c3 125 fLinkUpArray(0),
c978a69b 126 fTPCtracks(),
861b23c3 127 fEmcFirst(0),
128 fEmcLast(0),
129 fModule(0),
130 fTrackSegmentsInRun(0)
131{
132 // cpy ctor: no implementation yet
133 // requested by the Coding Convention
134 Fatal("cpy ctor", "not implemented") ;
135}
136
137
138//____________________________________________________________________________
139 AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
140{
141 // dtor
142 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
143 if (!fDefaultInit)
144 delete fLinkUpArray ;
861b23c3 145}
146
147//____________________________________________________________________________
148const TString AliPHOSTrackSegmentMakerv2::BranchName() const
149{
150
151 return GetName() ;
152}
153
154//____________________________________________________________________________
155void AliPHOSTrackSegmentMakerv2::FillOneModule()
156{
157 // Finds first and last indexes between which
158 // clusters from one PHOS module are
159
160 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
161 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
162
163 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
164
165 //First EMC clusters
166 Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
167 for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
168 ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
169 fEmcLast ++) ;
170
171 //Now TPC tracks
4b598bb1 172 if(fESD){
173 //Do it ones, only first time
174 if(fModule==1){
175 Int_t nTracks = fESD->GetNumberOfTracks();
c978a69b 176
177 Int_t nPHOSmod = geom->GetNModules() ;
178 if(fTPCtracks[0].size()<(UInt_t)nTracks){
179 for(Int_t i=0; i<nPHOSmod; i++)
180 fTPCtracks[i].resize(nTracks) ;
181 }
182 for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
4b598bb1 183 TVector3 inPHOS ;
184
185 //In this particular case we use fixed vertex position at zero
186 Double_t vtx[3]={0.,0.,0.} ;
187 AliESDtrack *track;
188 Double_t xyz[3] ;
4b598bb1 189 Double_t rEMC = geom->GetIPtoCrystalSurface() ; //Use here ideal geometry
190 for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
191 track = fESD->GetTrack(iTrack);
c978a69b 192 for(Int_t iTestMod=1; iTestMod<=nPHOSmod; iTestMod++){
193 Double_t modAngle=270.+geom->GetPHOSAngle(iTestMod) ;
194 modAngle*=TMath::Pi()/180. ;
195 track->Rotate(modAngle);
196 if (!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz))
197 continue; //track coord on the cylinder of PHOS radius
198 if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
199 continue;
200 //Check if this track hits PHOS
201 inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
202 Int_t modNum ;
203 Double_t x,z ;
204 geom->ImpactOnEmc(vtx, inPHOS.Theta(), inPHOS.Phi(), modNum, z, x) ;
205 if(modNum==iTestMod){
206 //Mark this track as one belonging to module
207 TrackInPHOS_t &t = fTPCtracks[modNum-1][fNtpcTracks[modNum-1]++] ;
208 t.track = track ;
209 t.x = x ;
210 t.z = z ;
211 break ;
212 }
4b598bb1 213 }
214 }
861b23c3 215 }
216 }
4b598bb1 217
861b23c3 218}
219
220//____________________________________________________________________________
221void AliPHOSTrackSegmentMakerv2::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
222 AliESDtrack *track,
223 Float_t &dx, Float_t &dz) const
224{
225 // Calculates the distance between the EMC RecPoint and the CPV RecPoint
226 // Clusters are sorted in "rows" and "columns" of width 1 cm
227
228 const AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
c978a69b 229 TVector3 emcGlobal; // Global position of the EMC recpoint
9ee9f78d 230 geom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);
4b598bb1 231
c978a69b 232 //Calculate actual distance to the PHOS surface
233 Double_t modAngle=270.+geom->GetPHOSAngle(emcClu->GetPHOSMod()) ;
234 modAngle*=TMath::Pi()/180. ;
235 Double_t rEMC = emcGlobal.X()*TMath::Cos(modAngle)+emcGlobal.Y()*TMath::Sin(modAngle) ;
236 track->Rotate(modAngle);
861b23c3 237 Double_t xyz[3] ;
c978a69b 238 if(!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)){
239 dx=999. ;
240 dz=999. ;
241 return ; //track coord on the cylinder of PHOS radius
861b23c3 242 }
c978a69b 243 if((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0){
861b23c3 244 dx=999. ;
245 dz=999. ;
c978a69b 246 return ;
247 }
248
249 dx=(emcGlobal.X()-xyz[0])*TMath::Sin(modAngle)-(emcGlobal.Y()-xyz[1])*TMath::Cos(modAngle) ;
250 dx=TMath::Sign(dx,(Float_t)(emcGlobal.X()-xyz[0])) ; //set direction
251 dz=emcGlobal.Z()-xyz[2] ;
861b23c3 252
253 return ;
254}
255//____________________________________________________________________________
256void AliPHOSTrackSegmentMakerv2::Init()
257{
258 // Make all memory allocations that are not possible in default constructor
259
260 AliPHOSGetter* gime = AliPHOSGetter::Instance();
261 if(!gime)
262 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
263
264 fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
265 if ( !gime->TrackSegmentMaker() ) {
266 gime->PostTrackSegmentMaker(this);
267 }
861b23c3 268}
269
270//____________________________________________________________________________
271void AliPHOSTrackSegmentMakerv2::InitParameters()
272{
273 //Initializes parameters
274 fRtpc = 4. ;
275 fEmcFirst = 0 ;
276 fEmcLast = 0 ;
277 fLinkUpArray = 0 ;
278 fWrite = kTRUE ;
279 fTrackSegmentsInRun = 0 ;
280 SetEventRange(0,-1) ;
281}
282
283
284//____________________________________________________________________________
c978a69b 285void AliPHOSTrackSegmentMakerv2::MakeLinks()
861b23c3 286{
287 // Finds distances (links) between all EMC and CPV clusters,
288 // which are not further apart from each other than fRcpv
289 // and sort them in accordance with this distance
290
291 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
292 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
293
294 fLinkUpArray->Clear() ;
295
296 AliPHOSEmcRecPoint * emcclu ;
297
298 Int_t iLinkUp = 0 ;
299
300 Int_t iEmcRP;
301 for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
302 emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
c978a69b 303 TVector3 vecEmc ;
304 emcclu->GetLocalPosition(vecEmc) ;
861b23c3 305 Int_t mod=emcclu->GetPHOSMod() ;
c978a69b 306 for(Int_t itr=0; itr<fNtpcTracks[mod-1] ; itr++){
307 TrackInPHOS_t &t = fTPCtracks[mod-1][itr] ;
308 //calculate raw distance
309 Double_t rawdx = t.x - vecEmc.X() ;
310 Double_t rawdz = t.z - vecEmc.Z() ;
311 if(TMath::Sqrt(rawdx*rawdx+rawdz*rawdz)<fRtpc){
312 Float_t dx,dz ;
313 //calculate presize difference accounting misalignment
314 GetDistanceInPHOSPlane(emcclu, t.track, dx,dz) ;
315 Int_t itrack = t.track->GetID() ;
861b23c3 316 new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(dx, dz, iEmcRP, itrack, -1) ;
317 }
318 }
319 }
861b23c3 320 fLinkUpArray->Sort() ; //first links with smallest distances
321}
322
323//____________________________________________________________________________
324void AliPHOSTrackSegmentMakerv2::MakePairs()
325{
326 // Using the previously made list of "links", we found the smallest link - i.e.
327 // link with the least distance between EMC and CPV and pointing to still
328 // unassigned RecParticles. We assign these RecPoints to TrackSegment and
329 // remove them from the list of "unassigned".
330
331 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
332
333 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
334 TClonesArray * trackSegments = gime->TrackSegments();
335
336 //Make arrays to mark clusters already chosen
337 Int_t * emcExist = 0;
338 if(fEmcLast > fEmcFirst)
339 emcExist = new Int_t[fEmcLast-fEmcFirst] ;
340
341 Int_t index;
342 for(index = 0; index <fEmcLast-fEmcFirst; index ++)
343 emcExist[index] = 1 ;
344
345
c978a69b 346 Bool_t * tpcExist = 0;
347 Int_t nTracks = fESD->GetNumberOfTracks();
861b23c3 348 if(nTracks>0)
c978a69b 349 tpcExist = new Bool_t[nTracks] ;
861b23c3 350
351 for(index = 0; index <nTracks; index ++)
352 tpcExist[index] = 1 ;
353
354
355 // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
356 TIter nextUp(fLinkUpArray) ;
357
358 AliPHOSLink * linkUp ;
359
360 AliPHOSCpvRecPoint * nullpointer = 0 ;
361
362 while ( (linkUp = static_cast<AliPHOSLink *>(nextUp()) ) ){
363
364 if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
365
366 if(tpcExist[linkUp->GetCpv()]){ //Track still exist
367 Float_t dx,dz ;
368 linkUp->GetXZ(dx,dz) ;
369 new ((* trackSegments)[fNTrackSegments])
370 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) ,
371 nullpointer,
372 linkUp->GetTrack(),dx,dz) ;
861b23c3 373 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
374 fNTrackSegments++ ;
375 emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found
376 //mark track as already used
377 tpcExist[linkUp->GetCpv()] = kFALSE ;
378 } //if CpvUp still exist
379 }
380 }
381
382 //look through emc recPoints left without CPV
383 if(emcExist){ //if there is emc rec point
384 Int_t iEmcRP ;
385 for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){
386 if(emcExist[iEmcRP] > 0 ){
387 new ((*trackSegments)[fNTrackSegments])
388 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP+fEmcFirst)),
389 nullpointer) ;
390 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
391 fNTrackSegments++;
392 }
393 }
394 }
395 delete [] emcExist ;
396 delete [] tpcExist ;
397}
398
399//____________________________________________________________________________
400void AliPHOSTrackSegmentMakerv2::Exec(Option_t *option)
401{
402 // Steering method to perform track segment construction for events
403 // in the range from fFirstEvent to fLastEvent.
404 // This range is optionally set by SetEventRange().
405 // if fLastEvent=-1 (by default), then process events until the end.
406
407 if(strstr(option,"tim"))
408 gBenchmark->Start("PHOSTSMaker");
409
410 if(strstr(option,"print")) {
411 Print() ;
412 return ;
413 }
414
415 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
416
417 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
418
419 if (fLastEvent == -1)
420 fLastEvent = gime->MaxEvent() - 1 ;
421 else
422 fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
423 Int_t nEvents = fLastEvent - fFirstEvent + 1;
424
425 Int_t ievent ;
426 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
427 gime->Event(ievent,"DR") ;
428 //Make some initializations
429 fNTrackSegments = 0 ;
430 fEmcFirst = 0 ;
431 fEmcLast = 0 ;
432
433 gime->TrackSegments()->Clear();
434
861b23c3 435 // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
436
437 for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
438 FillOneModule() ;
439 MakeLinks() ;
440 MakePairs() ;
441 }
442
443 WriteTrackSegments() ;
444
445 if(strstr(option,"deb"))
446 PrintTrackSegments(option);
447
448 //increment the total number of track segments per run
449 fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
450 }
451
452 if(strstr(option,"tim")){
453 gBenchmark->Stop("PHOSTSMaker");
454 Info("Exec", "took %f seconds for making TS %f seconds per event",
455 gBenchmark->GetCpuTime("PHOSTSMaker"),
456 gBenchmark->GetCpuTime("PHOSTSMaker")/nEvents) ;
457 }
458 if(fWrite) //do not unload in "on flight" mode
459 Unload();
460}
461//____________________________________________________________________________
861b23c3 462void AliPHOSTrackSegmentMakerv2::Unload()
463{
464 // Unloads the task from the folder
465 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
466 gime->PhosLoader()->UnloadRecPoints() ;
467 gime->PhosLoader()->UnloadTracks() ;
468}
469
470//____________________________________________________________________________
471void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
472{
473 // Print TrackSegmentMaker parameters
474
475 TString message("") ;
476 if( strcmp(GetName(), "") != 0 ) {
477 message = "\n======== AliPHOSTrackSegmentMakerv2 ========\n" ;
478 message += "Making Track segments\n" ;
479 message += "with parameters:\n" ;
480 message += " Maximal EMC - TPC distance (cm) %f\n" ;
481 message += "============================================\n" ;
482 Info("Print", message.Data(),fRtpc) ;
483 }
484 else
485 Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;
486}
487
488//____________________________________________________________________________
489void AliPHOSTrackSegmentMakerv2::WriteTrackSegments()
490{
491 // Writes found TrackSegments to TreeR. Creates branches
492 // "PHOSTS" and "AliPHOSTrackSegmentMaker" with the same title.
493 // In the former branch found TrackSegments are stored, while
494 // in the latter all parameters, with which TS were made.
495 // ROOT does not allow overwriting existing branches, therefore
496 // first we check, if branches with the same title already exist.
497 // If yes - exits without writing.
498
499 AliPHOSGetter *gime = AliPHOSGetter::Instance() ;
500
501 TClonesArray * trackSegments = gime->TrackSegments() ;
502 trackSegments->Expand(trackSegments->GetEntriesFast()) ;
503
504 if(fWrite){ //We write TreeT
505 TTree * treeT = gime->TreeT();
506
507 //First TS
508 Int_t bufferSize = 32000 ;
509 TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
510 tsBranch->Fill() ;
511
512 gime->WriteTracks("OVERWRITE");
513 gime->WriteTrackSegmentMaker("OVERWRITE");
514 }
515}
516
517
518//____________________________________________________________________________
519void AliPHOSTrackSegmentMakerv2::PrintTrackSegments(Option_t * option)
520{
521 // option deb - prints # of found TrackSegments
522 // option deb all - prints as well indexed of found RecParticles assigned to the TS
523
524 TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ;
525
526 Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
527 printf("nevent: %d\n", gAlice->GetEvNumber()) ;
528 printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
529
530 if(strstr(option,"all")) { // printing found TS
531 printf("TrackSegment # EMC RP# CPV RP#\n") ;
532 Int_t index;
533 for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
534 AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ;
535 printf(" %d %d %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetTrackIndex() ) ;
536 }
537 }
538}