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