e49794a0b742d6c2b44ddcd86bb03b52e99e7378
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv2.cxx
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$
20  *
21  */
22
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 
35 //   new TrackSegment. 
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.
40 // Use  case:
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") 
48 //                 
49 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) 
50 //
51
52 // --- ROOT system ---
53 #include "TTree.h"
54 #include "TBenchmark.h"
55
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"
64 #include "AliESD.h"
65 #include "AliESDtrack.h"
66
67 ClassImp( AliPHOSTrackSegmentMakerv2) 
68
69
70 //____________________________________________________________________________
71 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2() :
72   AliPHOSTrackSegmentMaker(),
73   fDefaultInit(kTRUE),
74   fWrite(kFALSE),
75   fNTrackSegments(0),
76   fRtpc(0.f),
77   fLinkUpArray(0),
78   fEmcFirst(0),
79   fEmcLast(0),
80   fModule(0),
81   fTrackSegmentsInRun(0)
82   
83 {
84   // default ctor (to be used mainly by Streamer)
85   InitParameters() ; 
86   for(Int_t i =0 ; i<5; i++)fTPC[i]=0 ;
87 }
88
89 //____________________________________________________________________________
90 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const TString & alirunFileName, const TString & eventFolderName) :
91   AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName),
92   fDefaultInit(kFALSE),
93   fWrite(kFALSE),
94   fNTrackSegments(0),
95   fRtpc(0.f),
96   fLinkUpArray(0),
97   fEmcFirst(0),
98   fEmcLast(0),
99   fModule(0),
100   fTrackSegmentsInRun(0)
101 {
102   // ctor
103   InitParameters() ; 
104   Init() ;
105   fESD = 0;
106 }
107
108 //____________________________________________________________________________
109 AliPHOSTrackSegmentMakerv2::AliPHOSTrackSegmentMakerv2(const AliPHOSTrackSegmentMakerv2 & tsm) :
110   AliPHOSTrackSegmentMaker(tsm),
111   fDefaultInit(kFALSE),
112   fWrite(kFALSE),
113   fNTrackSegments(0),
114   fRtpc(0.f),
115   fLinkUpArray(0),
116   fEmcFirst(0),
117   fEmcLast(0),
118   fModule(0),
119   fTrackSegmentsInRun(0)
120 {
121   // cpy ctor: no implementation yet
122   // requested by the Coding Convention
123   Fatal("cpy ctor", "not implemented") ;
124 }
125
126
127 //____________________________________________________________________________
128  AliPHOSTrackSegmentMakerv2::~AliPHOSTrackSegmentMakerv2()
129
130   // dtor
131   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
132   if (!fDefaultInit)  
133     delete fLinkUpArray ;
134   for(Int_t imod=0; imod<5; imod++){
135     if(fTPC[imod]) fTPC[imod]->Delete() ;
136   }
137 }
138
139 //____________________________________________________________________________
140 const TString AliPHOSTrackSegmentMakerv2::BranchName() const 
141 {  
142  
143   return GetName() ;
144 }
145
146 //____________________________________________________________________________
147 void  AliPHOSTrackSegmentMakerv2::FillOneModule()
148 {
149   // Finds first and last indexes between which 
150   // clusters from one PHOS module are
151
152   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
153   const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
154   
155   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
156  
157   //First EMC clusters
158   Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
159   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
160         ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
161       fEmcLast ++)  ;
162   
163   //Now TPC tracks
164   Int_t nTracks = fESD->GetNumberOfTracks();
165   TVector3 inPHOS ;
166  
167   //In this particular case we use fixed vertex position at zero
168   Double_t vtx[3]={0.,0.,0.} ;
169   AliESDtrack *track;
170   Double_t xyz[3] ;
171   Int_t nPHOSmod = geom->GetNModules() ;
172   for(Int_t imod=0 ; imod< nPHOSmod; imod++){
173     fTPC[imod]->Clear() ;
174   }
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)
181         continue;
182     //Check if this track hits PHOS
183     inPHOS.SetXYZ(xyz[0],xyz[1],xyz[2]);
184     Int_t modNum ; 
185     Double_t x,z ;
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) ;
190     }
191   } 
192       
193 }
194
195 //____________________________________________________________________________
196 void  AliPHOSTrackSegmentMakerv2::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,
197                                                          AliESDtrack *track,
198                                                          Float_t &dx, Float_t &dz) const
199 {
200   // Calculates the distance between the EMC RecPoint and the CPV RecPoint
201   // Clusters are sorted in "rows" and "columns" of width 1 cm
202
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 
207   Double_t xyz[3] ;
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] ;
212
213   }
214   else{
215     dx=999. ; 
216     dz=999. ;
217   }
218
219   return ;
220 }
221 //____________________________________________________________________________
222 void  AliPHOSTrackSegmentMakerv2::Init()
223 {
224   // Make all memory allocations that are not possible in default constructor
225   
226   AliPHOSGetter* gime = AliPHOSGetter::Instance();
227   if(!gime)
228     gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
229   
230   fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
231   if ( !gime->TrackSegmentMaker() ) {
232     gime->PostTrackSegmentMaker(this);
233   }
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() ;
238   }
239 }
240
241 //____________________________________________________________________________
242 void  AliPHOSTrackSegmentMakerv2::InitParameters()
243 {
244   //Initializes parameters
245   fRtpc      = 4. ;
246   fEmcFirst  = 0 ;    
247   fEmcLast   = 0 ;   
248   fLinkUpArray = 0 ;
249   fWrite                   = kTRUE ;
250   fTrackSegmentsInRun       = 0 ; 
251   SetEventRange(0,-1) ;
252 }
253
254
255 //____________________________________________________________________________
256 void  AliPHOSTrackSegmentMakerv2::MakeLinks()const
257
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
261   
262   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
263   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
264
265   fLinkUpArray->Clear() ;    
266
267   AliPHOSEmcRecPoint * emcclu ;
268
269   Int_t iLinkUp  = 0 ;
270   
271   Int_t iEmcRP;
272   for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
273     emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
274
275     Int_t mod=emcclu->GetPHOSMod() ;
276     TIter next(fTPC[mod-1]) ; 
277     AliESDtrack *track ;
278     Int_t itrack=0 ;
279     while((track= static_cast<AliESDtrack *>(next()))){
280       Float_t dx,dz ;
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) ;
284       }      
285     }
286   } 
287   
288   fLinkUpArray->Sort() ;  //first links with smallest distances
289 }
290
291 //____________________________________________________________________________
292 void  AliPHOSTrackSegmentMakerv2::MakePairs()
293
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". 
298
299   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
300
301   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
302   TClonesArray * trackSegments = gime->TrackSegments();
303     
304   //Make arrays to mark clusters already chosen
305   Int_t * emcExist = 0;
306   if(fEmcLast > fEmcFirst)
307     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
308   
309   Int_t index;
310   for(index = 0; index <fEmcLast-fEmcFirst; index ++)
311     emcExist[index] = 1 ;
312   
313   
314   Int_t * tpcExist = 0;
315   Int_t nTracks = fTPC[fModule]->GetSize() ;
316   if(nTracks>0)
317     tpcExist = new Int_t[nTracks] ;
318   
319   for(index = 0; index <nTracks; index ++)
320     tpcExist[index] = 1 ;
321   
322   
323   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
324   TIter nextUp(fLinkUpArray) ;
325   
326   AliPHOSLink * linkUp ;
327   
328   AliPHOSCpvRecPoint * nullpointer = 0 ;
329   
330   while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){  
331
332     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
333
334       if(tpcExist[linkUp->GetCpv()]){ //Track still exist
335          Float_t dx,dz ; 
336          linkUp->GetXZ(dx,dz) ;
337          new ((* trackSegments)[fNTrackSegments]) 
338            AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(linkUp->GetEmc())) , 
339                                nullpointer,
340                                linkUp->GetTrack(),dx,dz) ;
341          
342        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
343        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
348     } 
349   }        
350
351   //look through emc recPoints left without CPV
352   if(emcExist){ //if there is emc rec point
353     Int_t iEmcRP ;
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)), 
358                            nullpointer) ;
359        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
360        fNTrackSegments++;    
361       } 
362     }
363   }
364   delete [] emcExist ; 
365   delete [] tpcExist ; 
366 }
367
368 //____________________________________________________________________________
369 void  AliPHOSTrackSegmentMakerv2::Exec(Option_t *option)
370 {
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.
375   
376   if(strstr(option,"tim"))
377     gBenchmark->Start("PHOSTSMaker");
378  
379   if(strstr(option,"print")) {
380     Print() ; 
381     return ; 
382   }
383   
384   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
385  
386   const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
387
388   if (fLastEvent == -1) 
389     fLastEvent = gime->MaxEvent() - 1 ;
390   else 
391     fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
392   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
393
394   Int_t ievent ; 
395   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
396     gime->Event(ievent,"DR") ;
397    //Make some initializations 
398     fNTrackSegments = 0 ;
399     fEmcFirst = 0 ;    
400     fEmcLast  = 0 ;   
401     
402     gime->TrackSegments()->Clear();
403
404     GetVertex() ;
405     EvalRecPoints() ;
406
407     //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
408     
409     for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
410       FillOneModule() ; 
411       MakeLinks() ;
412       MakePairs() ;
413     }
414
415     WriteTrackSegments() ;
416
417     if(strstr(option,"deb"))
418       PrintTrackSegments(option);
419     
420     //increment the total number of track segments per run 
421     fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ; 
422   }
423   
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) ;
429    }
430   if(fWrite) //do not unload in "on flight" mode
431     Unload();
432 }
433 //____________________________________________________________________________
434 void AliPHOSTrackSegmentMakerv2::GetVertex(void)
435 { //extract vertex either using ESD or generator
436   
437   //Try to extract vertex from data
438   if(fESD){
439     const AliESDVertex *esdVtx = fESD->GetVertex() ;
440     if(esdVtx){
441       fVtx.SetXYZ(esdVtx->GetXv(),esdVtx->GetYv(),esdVtx->GetZv()) ;
442       return ;
443     }
444   }
445   
446   AliWarning("Can not read vertex from data, use fixed \n") ;
447   fVtx.SetXYZ(0.,0.,0.) ;
448  
449 }
450 //____________________________________________________________________________
451 void AliPHOSTrackSegmentMakerv2::EvalRecPoints(void)
452 { //calculate parameters of RecPoints using vertex and writing them
453
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) ;
461   }
462   emcRecPoints->Sort() ;
463
464   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
465   if(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) ;
469     }
470     cpvRecPoints->Sort() ;
471   }
472
473   //write recaculated RecPoints
474   gime->WriteRecPoints("OVERWRITE");
475
476 }
477 //____________________________________________________________________________
478 void AliPHOSTrackSegmentMakerv2::Unload() 
479 {
480   // Unloads the task from the folder
481   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
482   gime->PhosLoader()->UnloadRecPoints() ;
483   gime->PhosLoader()->UnloadTracks() ;
484 }
485
486 //____________________________________________________________________________
487 void AliPHOSTrackSegmentMakerv2::Print(const Option_t *)const
488 {
489   //  Print TrackSegmentMaker parameters
490
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) ;
499   }
500   else
501     Info("Print", "AliPHOSTrackSegmentMakerv2 not initialized ") ;
502 }
503
504 //____________________________________________________________________________
505 void AliPHOSTrackSegmentMakerv2::WriteTrackSegments()
506 {
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.
514
515   AliPHOSGetter *gime = AliPHOSGetter::Instance() ; 
516
517   TClonesArray * trackSegments = gime->TrackSegments() ; 
518   trackSegments->Expand(trackSegments->GetEntriesFast()) ;
519
520   if(fWrite){ //We write TreeT
521     TTree * treeT = gime->TreeT();
522     
523     //First TS
524     Int_t bufferSize = 32000 ; 
525     TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
526     tsBranch->Fill() ;  
527     
528     gime->WriteTracks("OVERWRITE");
529     gime->WriteTrackSegmentMaker("OVERWRITE");
530   }
531 }
532
533
534 //____________________________________________________________________________
535 void AliPHOSTrackSegmentMakerv2::PrintTrackSegments(Option_t * option)
536 {
537   // option deb - prints # of found TrackSegments
538   // option deb all - prints as well indexed of found RecParticles assigned to the TS
539
540   TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ; 
541
542   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
543   printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
544   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 
545   
546   if(strstr(option,"all")) {  // printing found TS
547     printf("TrackSegment #  EMC RP#  CPV RP#\n") ; 
548     Int_t index;
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() ) ; 
552     }   
553   }
554 }