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.3 2007/04/25 19:39:42 kharlov
21 * Track extracpolation improved
23 * Revision 1.2 2007/04/01 19:16:52 kharlov
24 * D.P.: Produce EMCTrackSegments using TPC/ITS tracks (no CPV)
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
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.
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")
55 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
58 // --- ROOT system ---
61 #include "TBenchmark.h"
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"
72 #include "AliESDtrack.h"
74 ClassImp( AliPHOSTrackSegmentMakerv2)
77 //____________________________________________________________________________
78 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2() :
79 AliPHOSTrackSegmentMaker(),
90 fTrackSegmentsInRun(0)
92 // default ctor (to be used mainly by Streamer)
96 //____________________________________________________________________________
97 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const TString & alirunFileName, const TString & eventFolderName) :
98 AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName),
109 fTrackSegmentsInRun(0)
117 //____________________________________________________________________________
118 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
119 AliPHOSTrackSegmentMaker(tsm),
120 fDefaultInit(kFALSE),
130 fTrackSegmentsInRun(0)
132 // cpy ctor: no implementation yet
133 // requested by the Coding Convention
134 Fatal("cpy ctor", "not implemented") ;
138 //____________________________________________________________________________
139 AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
142 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
144 delete fLinkUpArray ;
147 //____________________________________________________________________________
148 const TString AliPHOSTrackSegmentMakerv2::BranchName() const
154 //____________________________________________________________________________
155 void AliPHOSTrackSegmentMakerv2::FillOneModule()
157 // Finds first and last indexes between which
158 // clusters from one PHOS module are
160 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
161 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
163 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
166 Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
167 for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
168 ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
173 //Do it ones, only first time
175 Int_t nTracks = fESD->GetNumberOfTracks();
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) ;
182 for(Int_t i=0; i<5; i++)fNtpcTracks[i]=0 ;
185 //In this particular case we use fixed vertex position at zero
186 Double_t vtx[3]={0.,0.,0.} ;
189 Double_t rEMC = geom->GetIPtoCrystalSurface() ; //Use here ideal geometry
190 for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
191 track = fESD->GetTrack(iTrack);
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)
200 //Check if this track hits PHOS
201 inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
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]++] ;
220 //____________________________________________________________________________
221 void AliPHOSTrackSegmentMakerv2::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
223 Float_t &dx, Float_t &dz) const
225 // Calculates the distance between the EMC RecPoint and the CPV RecPoint
226 // Clusters are sorted in "rows" and "columns" of width 1 cm
228 const AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
229 TVector3 emcGlobal; // Global position of the EMC recpoint
230 geom->GetGlobalPHOS((AliPHOSRecPoint*)emcClu,emcGlobal);
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);
238 if(!track->GetXYZAt(rEMC, fESD->GetMagneticField(), xyz)){
241 return ; //track coord on the cylinder of PHOS radius
243 if((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0){
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] ;
255 //____________________________________________________________________________
256 void AliPHOSTrackSegmentMakerv2::Init()
258 // Make all memory allocations that are not possible in default constructor
260 AliPHOSGetter* gime = AliPHOSGetter::Instance();
262 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
264 fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
265 if ( !gime->TrackSegmentMaker() ) {
266 gime->PostTrackSegmentMaker(this);
270 //____________________________________________________________________________
271 void AliPHOSTrackSegmentMakerv2::InitParameters()
273 //Initializes parameters
279 fTrackSegmentsInRun = 0 ;
280 SetEventRange(0,-1) ;
284 //____________________________________________________________________________
285 void AliPHOSTrackSegmentMakerv2::MakeLinks()
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
291 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
292 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
294 fLinkUpArray->Clear() ;
296 AliPHOSEmcRecPoint * emcclu ;
301 for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
302 emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
304 emcclu->GetLocalPosition(vecEmc) ;
305 Int_t mod=emcclu->GetPHOSMod() ;
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){
313 //calculate presize difference accounting misalignment
314 GetDistanceInPHOSPlane(emcclu, t.track, dx,dz) ;
315 Int_t itrack = t.track->GetID() ;
316 new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(dx, dz, iEmcRP, itrack, -1) ;
320 fLinkUpArray->Sort() ; //first links with smallest distances
323 //____________________________________________________________________________
324 void AliPHOSTrackSegmentMakerv2::MakePairs()
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".
331 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
333 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
334 TClonesArray * trackSegments = gime->TrackSegments();
336 //Make arrays to mark clusters already chosen
337 Int_t * emcExist = 0;
338 if(fEmcLast > fEmcFirst)
339 emcExist = new Int_t[fEmcLast-fEmcFirst] ;
342 for(index = 0; index <fEmcLast-fEmcFirst; index ++)
343 emcExist[index] = 1 ;
346 Bool_t * tpcExist = 0;
347 Int_t nTracks = fESD->GetNumberOfTracks();
349 tpcExist = new Bool_t[nTracks] ;
351 for(index = 0; index <nTracks; index ++)
352 tpcExist[index] = 1 ;
355 // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
356 TIter nextUp(fLinkUpArray) ;
358 AliPHOSLink * linkUp ;
360 AliPHOSCpvRecPoint * nullpointer = 0 ;
362 while ( (linkUp = static_cast<AliPHOSLink *>(nextUp()) ) ){
364 if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
366 if(tpcExist[linkUp->GetCpv()]){ //Track still exist
368 linkUp->GetXZ(dx,dz) ;
369 new ((* trackSegments)[fNTrackSegments])
370 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) ,
372 linkUp->GetTrack(),dx,dz) ;
373 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(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
382 //look through emc recPoints left without CPV
383 if(emcExist){ //if there is emc rec point
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)),
390 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
399 //____________________________________________________________________________
400 void AliPHOSTrackSegmentMakerv2::Exec(Option_t *option)
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.
407 if(strstr(option,"tim"))
408 gBenchmark->Start("PHOSTSMaker");
410 if(strstr(option,"print")) {
415 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
417 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
419 if (fLastEvent == -1)
420 fLastEvent = gime->MaxEvent() - 1 ;
422 fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
423 Int_t nEvents = fLastEvent - fFirstEvent + 1;
426 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
427 gime->Event(ievent,"DR") ;
428 //Make some initializations
429 fNTrackSegments = 0 ;
433 gime->TrackSegments()->Clear();
435 // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
437 for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
443 WriteTrackSegments() ;
445 if(strstr(option,"deb"))
446 PrintTrackSegments(option);
448 //increment the total number of track segments per run
449 fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
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) ;
458 if(fWrite) //do not unload in "on flight" mode
461 //____________________________________________________________________________
462 void AliPHOSTrackSegmentMakerv2::Unload()
464 // Unloads the task from the folder
465 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
466 gime->PhosLoader()->UnloadRecPoints() ;
467 gime->PhosLoader()->UnloadTracks() ;
470 //____________________________________________________________________________
471 void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
473 // Print TrackSegmentMaker parameters
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) ;
485 Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;
488 //____________________________________________________________________________
489 void AliPHOSTrackSegmentMakerv2::WriteTrackSegments()
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.
499 AliPHOSGetter *gime = AliPHOSGetter::Instance() ;
501 TClonesArray * trackSegments = gime->TrackSegments() ;
502 trackSegments->Expand(trackSegments->GetEntriesFast()) ;
504 if(fWrite){ //We write TreeT
505 TTree * treeT = gime->TreeT();
508 Int_t bufferSize = 32000 ;
509 TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
512 gime->WriteTracks("OVERWRITE");
513 gime->WriteTrackSegmentMaker("OVERWRITE");
518 //____________________________________________________________________________
519 void AliPHOSTrackSegmentMakerv2::PrintTrackSegments(Option_t * option)
521 // option deb - prints # of found TrackSegments
522 // option deb all - prints as well indexed of found RecParticles assigned to the TS
524 TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ;
526 Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
527 printf("nevent: %d\n", gAlice->GetEvNumber()) ;
528 printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
530 if(strstr(option,"all")) { // printing found TS
531 printf("TrackSegment # EMC RP# CPV RP#\n") ;
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() ) ;