1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 /* History of cvs commits:
20 * Revision 1.73 2005/05/28 14:19:05 schutz
21 * Compilation warnings fixed by T.P.
25 //_________________________________________________________________________
26 // Implementation version 1 of algorithm class to construct PHOS track segments
27 // Track segment for PHOS is list of
28 // EMC RecPoint + (possibly) CPV RecPoint
29 // To find TrackSegments we do the following:
30 // for each EMC RecPoints we look at
31 // CPV RecPoints in the radious fRcpv.
32 // If there is such a CPV RecPoint,
33 // we make "Link" it is just indexes of EMC and CPV RecPoint and distance
34 // between them in the PHOS plane.
35 // Then we sort "Links" and starting from the
36 // least "Link" pointing to the unassined EMC and CPV RecPoints assing them to
38 // If there is no CPV RecPoint we make TrackSegment
39 // consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
40 //// In principle this class should be called from AliPHOSReconstructor, but
41 // one can use it as well in standalone mode.
43 // root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
44 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
45 // // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
46 // // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")
47 // root [1] t->ExecuteTask()
48 // root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
49 // root [4] t->ExecuteTask("deb all time")
51 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
54 // --- ROOT system ---
56 #include "TBenchmark.h"
58 // --- Standard library ---
59 #include "Riostream.h"
60 // --- AliRoot header files ---
61 #include "AliPHOSGeometry.h"
62 #include "AliPHOSTrackSegmentMakerv1.h"
63 #include "AliPHOSTrackSegment.h"
64 #include "AliPHOSLink.h"
65 #include "AliPHOSGetter.h"
67 #include "AliESDtrack.h"
69 ClassImp( AliPHOSTrackSegmentMakerv1)
72 //____________________________________________________________________________
73 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker()
75 // default ctor (to be used mainly by Streamer)
78 fDefaultInit = kTRUE ;
81 //____________________________________________________________________________
82 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
83 :AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName)
89 fDefaultInit = kFALSE ;
93 //____________________________________________________________________________
94 AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
97 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
103 //____________________________________________________________________________
104 const TString AliPHOSTrackSegmentMakerv1::BranchName() const
110 //____________________________________________________________________________
111 void AliPHOSTrackSegmentMakerv1::FillOneModule()
113 // Finds first and last indexes between which
114 // clusters from one PHOS module are
116 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
118 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
119 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
122 Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
123 for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
124 ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
128 Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
130 for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) &&
131 ((dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule );
136 //____________________________________________________________________________
137 Float_t AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSCpvRecPoint * cpvClu, Int_t &trackindex) const
139 // Calculates the distance between the EMC RecPoint and the CPV RecPoint
140 // Clusters are sorted in "rows" and "columns" of width 1 cm
142 //Float_t delta = 1 ; // Width of the rows in sorting of RecPoints (in cm)
143 // if you change this value, change it as well in xxxRecPoint::Compare()
144 Float_t distance2Cpv = fRcpv ;
145 Float_t distance2Track = fRtpc ;
147 trackindex = -1 ; // closest track within fRCpv
149 TVector3 vecEmc ; // Local position of EMC recpoint
150 TVector3 vecCpv ; // Local position of CPV recpoint propagated to EMC
151 TVector3 vecDist ; // Distance between local positions of two points
153 emcClu->GetLocalPosition(vecEmc) ;
154 cpvClu->GetLocalPosition(vecCpv) ;
157 if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){
159 // Find EMC-CPV distance
160 distance2Cpv = (vecCpv - vecEmc).Mag() ;
163 // Extrapolate the global track direction if any to CPV and find the closest track
164 Int_t nTracks = fESD->GetNumberOfTracks();
165 Int_t iClosestTrack = -1;
166 Double_t minDistance = 1e6;
167 Double_t pxyz[3], xyz[3];
169 AliPHOSGetter::Instance()->PHOSGeometry()->GetIPtoCrystalSurface();
171 for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
172 track = fESD->GetTrack(iTrack);
175 if (!track->GetXYZAt(rPHOS, fESD->GetMagneticField(), xyz))
176 continue; //track coord on the cylinder of PHOS radius
177 if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
179 if (!track->GetPxPyPzAt(rPHOS, fESD->GetMagneticField(), pxyz))
180 continue; // track momentum ibid.
181 vecDist = PropagateToPlane(xyz,pxyz,"CPV",cpvClu->GetPHOSMod());
182 // Info("GetDistanceInPHOSPlane","Track %d propagation to CPV = (%f,%f,%f)",
183 // iTrack,vecDist.X(),vecDist.Y(),vecDist.Z());
185 distance2Track = TMath::Sqrt(vecDist.X()*vecDist.X() + vecDist.Z()*vecDist.Z());
186 // Find the closest track to the EMC recpoint
187 if (distance2Track < minDistance) {
188 minDistance = distance2Track;
189 iClosestTrack = iTrack;
193 if (iClosestTrack != -1) {
194 track = fESD->GetTrack(iClosestTrack);
195 if (track->GetPxPyPzAt(rPHOS, fESD->GetMagneticField(), pxyz)) { // track momentum ibid.
196 TVector3 vecCpvGlobal; // Global position of the CPV recpoint
197 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
198 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
199 geom->GetGlobal((AliRecPoint*)cpvClu,vecCpvGlobal);
200 for (Int_t ixyz=0; ixyz<3; ixyz++)
201 xyz[ixyz] = vecCpvGlobal[ixyz];
202 vecDist = PropagateToPlane(xyz,pxyz,"EMC",cpvClu->GetPHOSMod());
203 // Info("GetDistanceInPHOSPlane","Track %d propagation to EMC = (%f,%f,%f)",
204 // iClosestTrack,vecDist.X(),vecDist.Y(),vecDist.Z());
206 distance2Track = TMath::Sqrt(vecDist.X()*vecDist.X() + vecDist.Z()*vecDist.Z());
210 // // If no ESD exists, than simply find EMC-CPV distance
211 // distance = (vecCpv - vecEmc).Mag() ;
213 //if(distance2Track < fRcpv + 2*delta )
214 if(distance2Track < fRtpc )
215 trackindex = iClosestTrack ;
218 // Info("GetDistanceInPHOSPlane","cpv-emc distance is %f cm",
222 return distance2Cpv ;
225 //____________________________________________________________________________
226 TVector3 AliPHOSTrackSegmentMakerv1::PropagateToPlane(Double_t *x, Double_t *p,
227 const char *det, Int_t moduleNumber) const
229 // Propagate a straight-line track from the origin point x
230 // along the direction p to the CPV or EMC module moduleNumber
231 // Returns a local position of such a propagation
233 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
234 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
235 TVector3 moduleCenter = geom->GetModuleCenter(det,moduleNumber);
237 TVector3 direction(p);
239 // Info("PropagateToCPV","Center of the %s module %d is (%f,%f,%f)",
240 // det,moduleNumber,moduleCenter[0],moduleCenter[1],moduleCenter[2]);
242 Double_t time = (moduleCenter.Mag2() - vertex.Dot(moduleCenter)) /
243 (direction.Dot(moduleCenter));
244 TVector3 globalIntersection = vertex + direction*time;
245 return geom->Global2Local(globalIntersection,moduleNumber);
248 //____________________________________________________________________________
249 void AliPHOSTrackSegmentMakerv1::Init()
251 // Make all memory allocations that are not possible in default constructor
253 AliPHOSGetter* gime = AliPHOSGetter::Instance();
255 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
257 fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
258 if ( !gime->TrackSegmentMaker() ) {
259 gime->PostTrackSegmentMaker(this);
263 //____________________________________________________________________________
264 void AliPHOSTrackSegmentMakerv1::InitParameters()
266 //Initializes parameters
275 fTrackSegmentsInRun = 0 ;
276 SetEventRange(0,-1) ;
280 //____________________________________________________________________________
281 void AliPHOSTrackSegmentMakerv1::MakeLinks()const
283 // Finds distances (links) between all EMC and CPV clusters,
284 // which are not further apart from each other than fRcpv
285 // and sort them in accordance with this distance
287 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
288 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
289 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
291 fLinkUpArray->Clear() ;
293 AliPHOSCpvRecPoint * cpv ;
294 AliPHOSEmcRecPoint * emcclu ;
299 for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
300 emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
304 for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) {
306 cpv = dynamic_cast<AliPHOSCpvRecPoint *>(cpvRecPoints->At(iCpv)) ;
308 Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, track) ;
312 new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(r, iEmcRP, iCpv, track) ;
317 fLinkUpArray->Sort() ; //first links with smallest distances
320 //____________________________________________________________________________
321 void AliPHOSTrackSegmentMakerv1::MakePairs()
323 // Using the previously made list of "links", we found the smallest link - i.e.
324 // link with the least distance between EMC and CPV and pointing to still
325 // unassigned RecParticles. We assign these RecPoints to TrackSegment and
326 // remove them from the list of "unassigned".
328 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
330 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
331 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
332 TClonesArray * trackSegments = gime->TrackSegments();
334 //Make arrays to mark clusters already chosen
335 Int_t * emcExist = 0;
336 if(fEmcLast > fEmcFirst)
337 emcExist = new Int_t[fEmcLast-fEmcFirst] ;
340 for(index = 0; index <fEmcLast-fEmcFirst; index ++)
341 emcExist[index] = 1 ;
343 Bool_t * cpvExist = 0;
344 if(fCpvLast > fCpvFirst)
345 cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
346 for(index = 0; index <fCpvLast-fCpvFirst; index ++)
347 cpvExist[index] = kTRUE ;
350 // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
351 TIter nextUp(fLinkUpArray) ;
353 AliPHOSLink * linkUp ;
355 AliPHOSCpvRecPoint * nullpointer = 0 ;
357 while ( (linkUp = static_cast<AliPHOSLink *>(nextUp()) ) ){
359 if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
361 if(cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist
362 new ((* trackSegments)[fNTrackSegments])
363 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) ,
364 dynamic_cast<AliPHOSCpvRecPoint *>(cpvRecPoints->At(linkUp->GetCpv())) ,
365 linkUp->GetTrack()) ;
367 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
369 emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found
370 //mark CPV recpoint as already used
371 cpvExist[linkUp->GetCpv()-fCpvFirst] = kFALSE ;
372 } //if CpvUp still exist
376 //look through emc recPoints left without CPV
377 if(emcExist){ //if there is emc rec point
379 for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){
380 if(emcExist[iEmcRP] > 0 ){
381 new ((*trackSegments)[fNTrackSegments])
382 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP+fEmcFirst)),
384 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
393 //____________________________________________________________________________
394 void AliPHOSTrackSegmentMakerv1::Exec(Option_t *option)
396 // Steering method to perform track segment construction for events
397 // in the range from fFirstEvent to fLastEvent.
398 // This range is optionally set by SetEventRange().
399 // if fLastEvent=-1 (by default), then process events until the end.
401 if(strstr(option,"tim"))
402 gBenchmark->Start("PHOSTSMaker");
404 if(strstr(option,"print")) {
409 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
411 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
413 if (fLastEvent == -1)
414 fLastEvent = gime->MaxEvent() - 1 ;
416 fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
417 Int_t nEvents = fLastEvent - fFirstEvent + 1;
420 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
421 gime->Event(ievent,"R") ;
422 //Make some initializations
423 fNTrackSegments = 0 ;
429 gime->TrackSegments()->Clear();
431 // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
433 for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
439 WriteTrackSegments() ;
441 if(strstr(option,"deb"))
442 PrintTrackSegments(option);
444 //increment the total number of track segments per run
445 fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
448 if(strstr(option,"tim")){
449 gBenchmark->Stop("PHOSTSMaker");
450 Info("Exec", "took %f seconds for making TS %f seconds per event",
451 gBenchmark->GetCpuTime("PHOSTSMaker"),
452 gBenchmark->GetCpuTime("PHOSTSMaker")/nEvents) ;
454 if(fWrite) //do not unload in "on flight" mode
458 //____________________________________________________________________________
459 void AliPHOSTrackSegmentMakerv1::Unload()
461 // Unloads the task from the folder
462 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
463 gime->PhosLoader()->UnloadRecPoints() ;
464 gime->PhosLoader()->UnloadTracks() ;
467 //____________________________________________________________________________
468 void AliPHOSTrackSegmentMakerv1::Print(const Option_t *)const
470 // Print TrackSegmentMaker parameters
472 TString message("") ;
473 if( strcmp(GetName(), "") != 0 ) {
474 message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ;
475 message += "Making Track segments\n" ;
476 message += "with parameters:\n" ;
477 message += " Maximal EMC - CPV distance (cm) %f\n" ;
478 message += "============================================\n" ;
479 Info("Print", message.Data(),fRcpv) ;
482 Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
485 //____________________________________________________________________________
486 void AliPHOSTrackSegmentMakerv1::WriteTrackSegments()
488 // Writes found TrackSegments to TreeR. Creates branches
489 // "PHOSTS" and "AliPHOSTrackSegmentMaker" with the same title.
490 // In the former branch found TrackSegments are stored, while
491 // in the latter all parameters, with which TS were made.
492 // ROOT does not allow overwriting existing branches, therefore
493 // first we check, if branches with the same title already exist.
494 // If yes - exits without writing.
496 AliPHOSGetter *gime = AliPHOSGetter::Instance() ;
498 TClonesArray * trackSegments = gime->TrackSegments() ;
499 trackSegments->Expand(trackSegments->GetEntriesFast()) ;
501 if(fWrite){ //We write TreeT
502 TTree * treeT = gime->TreeT();
505 Int_t bufferSize = 32000 ;
506 TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
509 gime->WriteTracks("OVERWRITE");
510 gime->WriteTrackSegmentMaker("OVERWRITE");
515 //____________________________________________________________________________
516 void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
518 // option deb - prints # of found TrackSegments
519 // option deb all - prints as well indexed of found RecParticles assigned to the TS
521 TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ;
523 Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
524 printf("nevent: %d\n", gAlice->GetEvNumber()) ;
525 printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
527 if(strstr(option,"all")) { // printing found TS
528 printf("TrackSegment # EMC RP# CPV RP#\n") ;
530 for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
531 AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ;
532 printf(" %d %d %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ;