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:
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
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.
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")
49 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH)
52 // --- ROOT system ---
54 #include "TBenchmark.h"
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"
65 #include "AliESDtrack.h"
67 ClassImp( AliPHOSTrackSegmentMakerv2)
70 //____________________________________________________________________________
71 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2() :
72 AliPHOSTrackSegmentMaker(),
81 fTrackSegmentsInRun(0)
84 // default ctor (to be used mainly by Streamer)
86 for(Int_t i =0 ; i<5; i++)fTPC[i]=0 ;
89 //____________________________________________________________________________
90 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const TString & alirunFileName, const TString & eventFolderName) :
91 AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName),
100 fTrackSegmentsInRun(0)
108 //____________________________________________________________________________
109 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
110 AliPHOSTrackSegmentMaker(tsm),
111 fDefaultInit(kFALSE),
119 fTrackSegmentsInRun(0)
121 // cpy ctor: no implementation yet
122 // requested by the Coding Convention
123 Fatal("cpy ctor", "not implemented") ;
127 //____________________________________________________________________________
128 AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
131 // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
133 delete fLinkUpArray ;
134 for(Int_t imod=0; imod<5; imod++){
135 if(fTPC[imod]) fTPC[imod]->Delete() ;
139 //____________________________________________________________________________
140 const TString AliPHOSTrackSegmentMakerv2::BranchName() const
146 //____________________________________________________________________________
147 void AliPHOSTrackSegmentMakerv2::FillOneModule()
149 // Finds first and last indexes between which
150 // clusters from one PHOS module are
152 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
153 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
155 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
158 Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
159 for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&
160 ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule );
164 Int_t nTracks = fESD->GetNumberOfTracks();
167 //In this particular case we use fixed vertex position at zero
168 Double_t vtx[3]={0.,0.,0.} ;
171 Int_t nPHOSmod = geom->GetNModules() ;
172 for(Int_t imod=0 ; imod< nPHOSmod; imod++){
173 fTPC[imod]->Clear() ;
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)
182 //Check if this track hits PHOS
183 inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
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) ;
195 //____________________________________________________________________________
196 void AliPHOSTrackSegmentMakerv2::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
198 Float_t &dx, Float_t &dz) const
200 // Calculates the distance between the EMC RecPoint and the CPV RecPoint
201 // Clusters are sorted in "rows" and "columns" of width 1 cm
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
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] ;
221 //____________________________________________________________________________
222 void AliPHOSTrackSegmentMakerv2::Init()
224 // Make all memory allocations that are not possible in default constructor
226 AliPHOSGetter* gime = AliPHOSGetter::Instance();
228 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
230 fLinkUpArray = new TClonesArray("AliPHOSLink", 1000);
231 if ( !gime->TrackSegmentMaker() ) {
232 gime->PostTrackSegmentMaker(this);
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() ;
241 //____________________________________________________________________________
242 void AliPHOSTrackSegmentMakerv2::InitParameters()
244 //Initializes parameters
250 fTrackSegmentsInRun = 0 ;
251 SetEventRange(0,-1) ;
255 //____________________________________________________________________________
256 void AliPHOSTrackSegmentMakerv2::MakeLinks()const
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
262 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
263 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
265 fLinkUpArray->Clear() ;
267 AliPHOSEmcRecPoint * emcclu ;
272 for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
273 emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
275 Int_t mod=emcclu->GetPHOSMod() ;
276 TIter next(fTPC[mod-1]) ;
279 while((track= static_cast<AliESDtrack *>(next()))){
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) ;
288 fLinkUpArray->Sort() ; //first links with smallest distances
291 //____________________________________________________________________________
292 void AliPHOSTrackSegmentMakerv2::MakePairs()
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".
299 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
301 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
302 TClonesArray * trackSegments = gime->TrackSegments();
304 //Make arrays to mark clusters already chosen
305 Int_t * emcExist = 0;
306 if(fEmcLast > fEmcFirst)
307 emcExist = new Int_t[fEmcLast-fEmcFirst] ;
310 for(index = 0; index <fEmcLast-fEmcFirst; index ++)
311 emcExist[index] = 1 ;
314 Int_t * tpcExist = 0;
315 Int_t nTracks = fTPC[fModule]->GetSize() ;
317 tpcExist = new Int_t[nTracks] ;
319 for(index = 0; index <nTracks; index ++)
320 tpcExist[index] = 1 ;
323 // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance
324 TIter nextUp(fLinkUpArray) ;
326 AliPHOSLink * linkUp ;
328 AliPHOSCpvRecPoint * nullpointer = 0 ;
330 while ( (linkUp = static_cast<AliPHOSLink *>(nextUp()) ) ){
332 if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
334 if(tpcExist[linkUp->GetCpv()]){ //Track still exist
336 linkUp->GetXZ(dx,dz) ;
337 new ((* trackSegments)[fNTrackSegments])
338 AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) ,
340 linkUp->GetTrack(),dx,dz) ;
342 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(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
351 //look through emc recPoints left without CPV
352 if(emcExist){ //if there is emc rec point
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)),
359 (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
368 //____________________________________________________________________________
369 void AliPHOSTrackSegmentMakerv2::Exec(Option_t *option)
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.
376 if(strstr(option,"tim"))
377 gBenchmark->Start("PHOSTSMaker");
379 if(strstr(option,"print")) {
384 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
386 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
388 if (fLastEvent == -1)
389 fLastEvent = gime->MaxEvent() - 1 ;
391 fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
392 Int_t nEvents = fLastEvent - fFirstEvent + 1;
395 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
396 gime->Event(ievent,"DR") ;
397 //Make some initializations
398 fNTrackSegments = 0 ;
402 gime->TrackSegments()->Clear();
407 // if(!ReadRecPoints(ievent)) continue; //reads RecPoints for event ievent
409 for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
415 WriteTrackSegments() ;
417 if(strstr(option,"deb"))
418 PrintTrackSegments(option);
420 //increment the total number of track segments per run
421 fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ;
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) ;
430 if(fWrite) //do not unload in "on flight" mode
433 //____________________________________________________________________________
434 void AliPHOSTrackSegmentMakerv2::GetVertex(void)
435 { //extract vertex either using ESD or generator
437 //Try to extract vertex from data
439 const AliESDVertex *esdVtx = fESD->GetVertex() ;
441 fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
446 AliWarning("Can not read vertex from data, use fixed \n") ;
447 fVtx.SetXYZ(0.,0.,0.) ;
450 //____________________________________________________________________________
451 void AliPHOSTrackSegmentMakerv2::EvalRecPoints(void)
452 { //calculate parameters of RecPoints using vertex and writing them
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) ;
462 emcRecPoints->Sort() ;
464 TObjArray * cpvRecPoints = gime->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) ;
470 cpvRecPoints->Sort() ;
473 //write recaculated RecPoints
474 gime->WriteRecPoints("OVERWRITE");
477 //____________________________________________________________________________
478 void AliPHOSTrackSegmentMakerv2::Unload()
480 // Unloads the task from the folder
481 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
482 gime->PhosLoader()->UnloadRecPoints() ;
483 gime->PhosLoader()->UnloadTracks() ;
486 //____________________________________________________________________________
487 void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
489 // Print TrackSegmentMaker parameters
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) ;
501 Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;
504 //____________________________________________________________________________
505 void AliPHOSTrackSegmentMakerv2::WriteTrackSegments()
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.
515 AliPHOSGetter *gime = AliPHOSGetter::Instance() ;
517 TClonesArray * trackSegments = gime->TrackSegments() ;
518 trackSegments->Expand(trackSegments->GetEntriesFast()) ;
520 if(fWrite){ //We write TreeT
521 TTree * treeT = gime->TreeT();
524 Int_t bufferSize = 32000 ;
525 TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
528 gime->WriteTracks("OVERWRITE");
529 gime->WriteTrackSegmentMaker("OVERWRITE");
534 //____________________________________________________________________________
535 void AliPHOSTrackSegmentMakerv2::PrintTrackSegments(Option_t * option)
537 // option deb - prints # of found TrackSegments
538 // option deb all - prints as well indexed of found RecParticles assigned to the TS
540 TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ;
542 Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ;
543 printf("nevent: %d\n", gAlice->GetEvNumber()) ;
544 printf(" Found %d TrackSegments\n", trackSegments->GetEntriesFast() );
546 if(strstr(option,"all")) { // printing found TS
547 printf("TrackSegment # EMC RP# CPV RP#\n") ;
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() ) ;