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.85 2007/03/28 19:18:15 kharlov
21 * RecPoints recalculation in TSM removed
23 * Revision 1.84 2007/03/07 07:01:21 hristov
24 * Fixing copy/paste erro. Additional protections
26 * Revision 1.83 2007/03/06 21:07:37 kharlov
27 * DP: xz CPV-EMC distance filled to TS
29 * Revision 1.82 2007/03/06 06:54:48 kharlov
30 * DP:Calculation of cluster properties dep. on vertex added
32 * Revision 1.81 2007/02/05 10:02:40 kharlov
33 * Module numbering is corrected
35 * Revision 1.80 2006/08/28 10:01:56 kharlov
36 * Effective C++ warnings fixed (Timur Pocheptsov)
38 * Revision 1.79 2006/04/25 12:41:15 hristov
39 * Moving non-persistent data to AliESDfriend (Yu.Belikov)
41 * Revision 1.78 2005/11/18 13:04:51 hristov
44 * Revision 1.77 2005/11/17 23:34:36 hristov
47 * Revision 1.76 2005/11/17 22:29:12 hristov
48 * Faster version, no attempt to match tracks outside the PHOS acceptance
50 * Revision 1.75 2005/11/17 12:35:27 hristov
51 * Use references instead of objects. Avoid to create objects when they are not really needed
53 * Revision 1.74 2005/07/08 14:01:36 hristov
54 * Tracking in non-uniform nmagnetic field (Yu.Belikov)
56 * Revision 1.73 2005/05/28 14:19:05 schutz
57 * Compilation warnings fixed by T.P.
61 //_________________________________________________________________________
62 // Implementation version 1 of algorithm class to construct PHOS track segments
63 // Track segment for PHOS is list of
64 // EMC RecPoint + (possibly) CPV RecPoint
65 // To find TrackSegments we do the following:
66 // for each EMC RecPoints we look at
67 // CPV RecPoints in the radious fRcpv.
68 // If there is such a CPV RecPoint,
69 // we make "Link" it is just indexes of EMC and CPV RecPoint and distance
70 // between them in the PHOS plane.
71 // Then we sort "Links" and starting from the
72 // least "Link" pointing to the unassined EMC and CPV RecPoints assing them to
74 // If there is no CPV RecPoint we make TrackSegment
75 // consisting from EMC alone. There is no TrackSegments without EMC RecPoint.
76 //// In principle this class should be called from AliPHOSReconstructor, but
77 // one can use it as well in standalone mode.
79 // root [0] AliPHOSTrackSegmentMakerv1 * t = new AliPHOSTrackSegmentMaker("galice.root", "tracksegmentsname", "recpointsname")
80 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
81 // // reads gAlice from header file "galice.root", uses recpoints stored in the branch names "recpointsname" (default = "Default")
82 // // and saves recpoints in branch named "tracksegmentsname" (default = "recpointsname")
83 // root [1] t->ExecuteTask()
84 // root [3] t->SetTrackSegmentsBranch("max distance 5 cm")
85 // root [4] t->ExecuteTask("deb all time")
87 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
90 // --- ROOT system ---
92 #include "TBenchmark.h"
94 // --- Standard library ---
95 #include "Riostream.h"
96 // --- AliRoot header files ---
97 #include "AliPHOSGeometry.h"
98 #include "AliPHOSTrackSegmentMakerv1.h"
99 #include "AliPHOSTrackSegment.h"
100 #include "AliPHOSLink.h"
101 #include "AliPHOSGetter.h"
103 #include "AliESDtrack.h"
105 ClassImp( AliPHOSTrackSegmentMakerv1)
108 //____________________________________________________________________________
109 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() :
110 AliPHOSTrackSegmentMaker(),
122 fTrackSegmentsInRun(0)
125 // default ctor (to be used mainly by Streamer)
129 //____________________________________________________________________________
130 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const TString & alirunFileName, const TString & eventFolderName) :
131 AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName),
132 fDefaultInit(kFALSE),
143 fTrackSegmentsInRun(0)
152 AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const AliPHOSTrackSegmentMakerv1 & tsm) :
153 AliPHOSTrackSegmentMaker(tsm),
154 fDefaultInit(kFALSE),
165 fTrackSegmentsInRun(0)
167 // cpy ctor: no implementation yet
168 // requested by the Coding Convention
169 Fatal("cpy ctor", "not implemented") ;
173 //____________________________________________________________________________
174 AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
177 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
179 delete fLinkUpArray ;
183 //____________________________________________________________________________
184 const TString AliPHOSTrackSegmentMakerv1::BranchName() const
190 //____________________________________________________________________________
191 void AliPHOSTrackSegmentMakerv1::FillOneModule()
193 // Finds first and last indexes between which
194 // clusters from one PHOS module are
196 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
198 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
199 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
202 Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
203 for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
204 ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
208 Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
210 for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) &&
211 ((dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule );
216 //____________________________________________________________________________
217 void AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
218 AliPHOSCpvRecPoint * cpvClu,
220 Float_t &dx, Float_t &dz) const
222 // Calculates the distance between the EMC RecPoint and the CPV RecPoint
223 // Clusters are sorted in "rows" and "columns" of width 1 cm
225 // Float_t delta = 1 ; // Width of the rows in sorting of RecPoints (in cm)
226 // // if you change this value, change it as well in xxxRecPoint::Compare()
227 Float_t distance2Track = fRtpc ;
229 trackindex = -1 ; // closest track within fRCpv
231 TVector3 vecEmc ; // Local position of EMC recpoint
232 TVector3 vecP ; // Momentum direction at CPV plain
233 TVector3 vecPloc ; // Momentum direction at CPV plain
236 if(emcClu->GetPHOSMod() != cpvClu->GetPHOSMod()){
242 emcClu->GetLocalPosition(vecEmc) ;
244 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
245 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
247 Double_t xCPV,zCPV ; //EMC-projected coordinates of CPV cluster
248 TVector3 cpvGlobal; // Global position of the CPV recpoint
249 geom->GetGlobal((AliRecPoint*)cpvClu,cpvGlobal);
250 Double_t vtxCPV[3]={cpvGlobal.X(),cpvGlobal.Y(),cpvGlobal.Z()} ;
253 //if no track information available, assume straight line from IP to emcal
255 geom->ImpactOnEmc(vtxCPV,cpvGlobal.Theta(),cpvGlobal.Phi(),dummyMod,xCPV,zCPV) ;
256 dx=xCPV - vecEmc.X() ;
257 dz=zCPV - vecEmc.Z() ;
261 //if there is ESD try to correct distance using TPC information on particle direct in CPV
264 Double_t rCPV = cpvGlobal.Pt() ;// Radius from IP to current point
266 // Extrapolate the global track direction if any to CPV and find the closest track
267 Int_t nTracks = fESD->GetNumberOfTracks();
268 Int_t iClosestTrack = -1;
269 Double_t minDistance = 1.e6;
275 for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
276 track = fESD->GetTrack(iTrack);
277 if (!track->GetXYZAt(rCPV, fESD->GetMagneticField(), xyz))
278 continue; //track coord on the cylinder of PHOS radius
279 if ((TMath::Abs(xyz[0])+TMath::Abs(xyz[1])+TMath::Abs(xyz[2]))<=0)
281 //Check if this track hits PHOS
282 inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
283 distance2Track = inPHOS.Angle(cpvGlobal) ;
284 // Find the closest track to the CPV recpoint
285 if (distance2Track < minDistance) {
286 minDistance = distance2Track;
287 iClosestTrack = iTrack;
291 if (iClosestTrack != -1) {
292 track = fESD->GetTrack(iClosestTrack);
293 if (track->GetPxPyPzAt(rCPV, fESD->GetMagneticField(), pxyz)) { // track momentum ibid.
294 vecP.SetXYZ(pxyz[0],pxyz[1],pxyz[2]);
296 geom->ImpactOnEmc(vtxCPV,vecP.Theta(),vecP.Phi(),dummyMod,xCPV,zCPV) ;
300 if(minDistance < fRtpc ){
301 trackindex = iClosestTrack ;
305 dx=xCPV - vecEmc.X() ;
306 dz=zCPV - vecEmc.Z() ;
314 //____________________________________________________________________________
315 void AliPHOSTrackSegmentMakerv1::Init()
317 // Make all memory allocations that are not possible in default constructor
319 AliPHOSGetter* gime = AliPHOSGetter::Instance();
321 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
323 fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
324 if ( !gime->TrackSegmentMaker() ) {
325 gime->PostTrackSegmentMaker(this);
329 //____________________________________________________________________________
330 void AliPHOSTrackSegmentMakerv1::InitParameters()
332 //Initializes parameters
341 fTrackSegmentsInRun = 0 ;
342 SetEventRange(0,-1) ;
346 //____________________________________________________________________________
347 void AliPHOSTrackSegmentMakerv1::MakeLinks()const
349 // Finds distances (links) between all EMC and CPV clusters,
350 // which are not further apart from each other than fRcpv
351 // and sort them in accordance with this distance
353 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
354 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
355 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
357 fLinkUpArray->Clear() ;
359 AliPHOSCpvRecPoint * cpv ;
360 AliPHOSEmcRecPoint * emcclu ;
365 for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
366 emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
370 for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) {
372 cpv = dynamic_cast<AliPHOSCpvRecPoint *>(cpvRecPoints->At(iCpv)) ;
375 GetDistanceInPHOSPlane(emcclu, cpv, track,dx,dz) ;
376 if(TMath::Sqrt(dx*dx+dz*dz) < fRcpv ){
377 new ((*fLinkUpArray)[iLinkUp++]) AliPHOSLink(dx, dz, iEmcRP, iCpv, track) ;
382 fLinkUpArray->Sort() ; //first links with smallest distances
385 //____________________________________________________________________________
386 void AliPHOSTrackSegmentMakerv1::MakePairs()
388 // Using the previously made list of "links", we found the smallest link - i.e.
389 // link with the least distance between EMC and CPV and pointing to still
390 // unassigned RecParticles. We assign these RecPoints to TrackSegment and
391 // remove them from the list of "unassigned".
393 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
395 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
396 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
397 TClonesArray * trackSegments = gime->TrackSegments();
399 //Make arrays to mark clusters already chosen
400 Int_t * emcExist = 0;
401 if(fEmcLast > fEmcFirst)
402 emcExist = new Int_t[fEmcLast-fEmcFirst] ;
405 for(index = 0; index <fEmcLast-fEmcFirst; index ++)
406 emcExist[index] = 1 ;
408 Bool_t * cpvExist = 0;
409 if(fCpvLast > fCpvFirst)
410 cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
411 for(index = 0; index <fCpvLast-fCpvFirst; index ++)
412 cpvExist[index] = kTRUE ;
415 // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
416 TIter nextUp(fLinkUpArray) ;
418 AliPHOSLink * linkUp ;
420 AliPHOSCpvRecPoint * nullpointer = 0 ;
422 while ( (linkUp = static_cast<AliPHOSLink *>(nextUp()) ) ){
424 if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
426 if(cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist
428 linkUp->GetXZ(dx,dz) ;
429 new ((* trackSegments)[fNTrackSegments])
430 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) ,
431 dynamic_cast<AliPHOSCpvRecPoint *>(cpvRecPoints->At(linkUp->GetCpv())) ,
432 linkUp->GetTrack(),dx,dz) ;
434 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
436 emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc that Cpv was found
437 //mark CPV recpoint as already used
438 cpvExist[linkUp->GetCpv()-fCpvFirst] = kFALSE ;
439 } //if CpvUp still exist
443 //look through emc recPoints left without CPV
444 if(emcExist){ //if there is emc rec point
446 for(iEmcRP = 0; iEmcRP < fEmcLast-fEmcFirst ; iEmcRP++ ){
447 if(emcExist[iEmcRP] > 0 ){
448 new ((*trackSegments)[fNTrackSegments])
449 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP+fEmcFirst)),
451 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
460 //____________________________________________________________________________
461 void AliPHOSTrackSegmentMakerv1::Exec(Option_t *option)
463 // Steering method to perform track segment construction for events
464 // in the range from fFirstEvent to fLastEvent.
465 // This range is optionally set by SetEventRange().
466 // if fLastEvent=-1 (by default), then process events until the end.
468 if(strstr(option,"tim"))
469 gBenchmark->Start("PHOSTSMaker");
471 if(strstr(option,"print")) {
476 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
478 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
480 if (fLastEvent == -1)
481 fLastEvent = gime->MaxEvent() - 1 ;
483 fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
484 Int_t nEvents = fLastEvent - fFirstEvent + 1;
487 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
488 gime->Event(ievent,"DR") ;
489 //Make some initializations
490 fNTrackSegments = 0 ;
496 gime->TrackSegments()->Clear();
498 // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
500 for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
506 WriteTrackSegments() ;
508 if(strstr(option,"deb"))
509 PrintTrackSegments(option);
511 //increment the total number of track segments per run
512 fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
515 if(strstr(option,"tim")){
516 gBenchmark->Stop("PHOSTSMaker");
517 Info("Exec", "took %f seconds for making TS %f seconds per event",
518 gBenchmark->GetCpuTime("PHOSTSMaker"),
519 gBenchmark->GetCpuTime("PHOSTSMaker")/nEvents) ;
521 if(fWrite) //do not unload in "on flight" mode
524 //____________________________________________________________________________
525 void AliPHOSTrackSegmentMakerv1::Unload()
527 // Unloads the task from the folder
528 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
529 gime->PhosLoader()->UnloadRecPoints() ;
530 gime->PhosLoader()->UnloadTracks() ;
533 //____________________________________________________________________________
534 void AliPHOSTrackSegmentMakerv1::Print(const Option_t *)const
536 // Print TrackSegmentMaker parameters
538 TString message("") ;
539 if( strcmp(GetName(), "") != 0 ) {
540 message = "\n======== AliPHOSTrackSegmentMakerv1 ========\n" ;
541 message += "Making Track segments\n" ;
542 message += "with parameters:\n" ;
543 message += " Maximal EMC - CPV distance (cm) %f\n" ;
544 message += "============================================\n" ;
545 Info("Print", message.Data(),fRcpv) ;
548 Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
551 //____________________________________________________________________________
552 void AliPHOSTrackSegmentMakerv1::WriteTrackSegments()
554 // Writes found TrackSegments to TreeR. Creates branches
555 // "PHOSTS" and "AliPHOSTrackSegmentMaker" with the same title.
556 // In the former branch found TrackSegments are stored, while
557 // in the latter all parameters, with which TS were made.
558 // ROOT does not allow overwriting existing branches, therefore
559 // first we check, if branches with the same title already exist.
560 // If yes - exits without writing.
562 AliPHOSGetter *gime = AliPHOSGetter::Instance() ;
564 TClonesArray * trackSegments = gime->TrackSegments() ;
565 trackSegments->Expand(trackSegments->GetEntriesFast()) ;
567 if(fWrite){ //We write TreeT
568 TTree * treeT = gime->TreeT();
571 Int_t bufferSize = 32000 ;
572 TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
575 gime->WriteTracks("OVERWRITE");
576 gime->WriteTrackSegmentMaker("OVERWRITE");
581 //____________________________________________________________________________
582 void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
584 // option deb - prints # of found TrackSegments
585 // option deb all - prints as well indexed of found RecParticles assigned to the TS
587 TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ;
589 Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
590 printf("nevent: %d\n", AliPHOSGetter::Instance()->EventNumber()) ;
591 printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
593 if(strstr(option,"all")) { // printing found TS
594 printf("TrackSegment # EMC RP# CPV RP#\n") ;
596 for (index = 0 ; index <trackSegments->GetEntriesFast() ; index++) {
597 AliPHOSTrackSegment * ts = (AliPHOSTrackSegment * )trackSegments->At(index) ;
598 printf(" %d %d %d \n", ts->GetIndexInList(), ts->GetEmcIndex(), ts->GetCpvIndex() ) ;