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