]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PHOS/AliPHOSTrackSegmentMakerv1.cxx
Decalibrate simulated data
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv1.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  * Revision 1.73  2005/05/28 14:19:05  schutz
21  * Compilation warnings fixed by T.P.
22  *
23  */
24
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 
37 //   new TrackSegment. 
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.
42 // Use  case:
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") 
50 //                 
51 //*-- Author: Dmitri Peressounko (RRC Ki & SUBATECH) & Yves Schutz (SUBATECH) 
52 //
53
54 // --- ROOT system ---
55 #include "TTree.h"
56 #include "TBenchmark.h"
57
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"
66 #include "AliESD.h"
67 #include "AliESDtrack.h"
68
69 ClassImp( AliPHOSTrackSegmentMakerv1) 
70
71
72 //____________________________________________________________________________
73   AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() : AliPHOSTrackSegmentMaker()
74 {
75   // default ctor (to be used mainly by Streamer)
76
77   InitParameters() ; 
78   fDefaultInit = kTRUE ; 
79 }
80
81 //____________________________________________________________________________
82  AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const TString alirunFileName, const TString eventFolderName)
83    :AliPHOSTrackSegmentMaker(alirunFileName, eventFolderName)
84 {
85   // ctor
86
87   InitParameters() ; 
88   Init() ;
89   fDefaultInit = kFALSE ; 
90   fESD = 0;
91 }
92
93 //____________________________________________________________________________
94  AliPHOSTrackSegmentMakerv1::~AliPHOSTrackSegmentMakerv1()
95
96   // dtor
97   // fDefaultInit = kTRUE if TrackSegmentMaker created by default ctor (to get just the parameters)
98   if (!fDefaultInit)  
99     delete fLinkUpArray ;
100 }
101
102
103 //____________________________________________________________________________
104 const TString AliPHOSTrackSegmentMakerv1::BranchName() const 
105 {  
106  
107   return GetName() ;
108 }
109
110 //____________________________________________________________________________
111 void  AliPHOSTrackSegmentMakerv1::FillOneModule()
112 {
113   // Finds first and last indexes between which 
114   // clusters from one PHOS module are
115
116   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
117   
118   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
119   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
120  
121   //First EMC clusters
122   Int_t totalEmc = emcRecPoints->GetEntriesFast() ;
123   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
124         ((dynamic_cast<AliPHOSRecPoint *>(emcRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
125       fEmcLast ++)  ;
126   
127   //Now CPV clusters
128   Int_t totalCpv = cpvRecPoints->GetEntriesFast() ;
129
130     for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && 
131          ((dynamic_cast<AliPHOSRecPoint *>(cpvRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); 
132        fCpvLast ++) ;
133       
134 }
135
136 //____________________________________________________________________________
137 Float_t  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * emcClu,AliPHOSCpvRecPoint * cpvClu, Int_t &trackindex) const
138 {
139   // Calculates the distance between the EMC RecPoint and the CPV RecPoint
140   // Clusters are sorted in "rows" and "columns" of width 1 cm
141
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 ; 
146
147   trackindex = -1 ; // closest track within fRCpv 
148
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
152   
153   emcClu->GetLocalPosition(vecEmc) ;
154   cpvClu->GetLocalPosition(vecCpv) ;
155
156   //toofar = kTRUE ;
157   if(emcClu->GetPHOSMod() == cpvClu->GetPHOSMod()){ 
158
159     // Find EMC-CPV distance
160     distance2Cpv = (vecCpv - vecEmc).Mag() ;
161     
162     if (fESD != 0x0) {
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];
168       Double_t rPHOS = 
169         AliPHOSGetter::Instance()->PHOSGeometry()->GetIPtoCrystalSurface();
170       AliESDtrack *track;
171       for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
172         track = fESD->GetTrack(iTrack);
173         if (track->IsPHOS()) 
174           continue ; 
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)
178            continue;
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());
184         vecDist -= vecCpv;
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;
190         }
191       }
192
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());
205         vecDist -= vecEmc;
206         distance2Track = TMath::Sqrt(vecDist.X()*vecDist.X() + vecDist.Z()*vecDist.Z());
207         }
208       }
209 //     } else {
210 //       // If no ESD exists, than simply find EMC-CPV distance
211 //       distance = (vecCpv - vecEmc).Mag() ;
212     
213       //if(distance2Track < fRcpv + 2*delta )
214       if(distance2Track < fRtpc )
215         trackindex = iClosestTrack ; 
216       //      toofar = kFALSE ;
217     }
218     //     Info("GetDistanceInPHOSPlane","cpv-emc distance is %f cm",
219     //   distance);
220   }
221   
222   return distance2Cpv ;
223 }
224
225 //____________________________________________________________________________
226 TVector3  AliPHOSTrackSegmentMakerv1::PropagateToPlane(Double_t *x, Double_t *p,
227                                                        const char *det, Int_t moduleNumber) const
228 {
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
232
233   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
234   const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
235   TVector3 moduleCenter = geom->GetModuleCenter(det,moduleNumber);
236   TVector3 vertex(x);
237   TVector3 direction(p);
238
239 //   Info("PropagateToCPV","Center of the %s module %d is (%f,%f,%f)",
240 //        det,moduleNumber,moduleCenter[0],moduleCenter[1],moduleCenter[2]);
241
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);
246 }
247
248 //____________________________________________________________________________
249 void  AliPHOSTrackSegmentMakerv1::Init()
250 {
251   // Make all memory allocations that are not possible in default constructor
252   
253   AliPHOSGetter* gime = AliPHOSGetter::Instance();
254   if(!gime)
255     gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
256   
257   fLinkUpArray  = new TClonesArray("AliPHOSLink", 1000); 
258   if ( !gime->TrackSegmentMaker() ) {
259     gime->PostTrackSegmentMaker(this);
260   }
261 }
262
263 //____________________________________________________________________________
264 void  AliPHOSTrackSegmentMakerv1::InitParameters()
265 {
266   //Initializes parameters
267   fRcpv      = 10. ;
268   fRtpc      = 4. ;
269   fEmcFirst  = 0 ;    
270   fEmcLast   = 0 ;   
271   fCpvFirst  = 0 ;   
272   fCpvLast   = 0 ;   
273   fLinkUpArray = 0 ;
274   fWrite                   = kTRUE ;
275   fTrackSegmentsInRun       = 0 ; 
276   SetEventRange(0,-1) ;
277 }
278
279
280 //____________________________________________________________________________
281 void  AliPHOSTrackSegmentMakerv1::MakeLinks()const
282
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
286   
287   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
288   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
289   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
290
291   fLinkUpArray->Clear() ;    
292
293   AliPHOSCpvRecPoint * cpv ;
294   AliPHOSEmcRecPoint * emcclu ;
295
296   Int_t iLinkUp  = 0 ;
297   
298   Int_t iEmcRP;
299   for(iEmcRP = fEmcFirst; iEmcRP < fEmcLast; iEmcRP++ ) {
300     emcclu = dynamic_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(iEmcRP)) ;
301
302     //Bool_t toofar ;        
303     Int_t iCpv = 0 ;    
304     for(iCpv = fCpvFirst; iCpv < fCpvLast;iCpv++ ) { 
305       
306       cpv = dynamic_cast<AliPHOSCpvRecPoint *>(cpvRecPoints->At(iCpv)) ;
307       Int_t track = -1 ; 
308       Float_t r = GetDistanceInPHOSPlane(emcclu, cpv, track) ;     
309       //      if(toofar)
310       //        continue ;       
311       if(r < fRcpv) { 
312         new ((*fLinkUpArray)[iLinkUp++])  AliPHOSLink(r, iEmcRP, iCpv, track) ;
313       }      
314     }
315   } 
316   
317   fLinkUpArray->Sort() ;  //first links with smallest distances
318 }
319
320 //____________________________________________________________________________
321 void  AliPHOSTrackSegmentMakerv1::MakePairs()
322
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". 
327
328   AliPHOSGetter * gime = AliPHOSGetter::Instance() ; 
329
330   TObjArray * emcRecPoints = gime->EmcRecPoints() ; 
331   TObjArray * cpvRecPoints = gime->CpvRecPoints() ; 
332   TClonesArray * trackSegments = gime->TrackSegments();
333     
334   //Make arrays to mark clusters already chosen
335   Int_t * emcExist = 0;
336   if(fEmcLast > fEmcFirst)
337     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
338   
339   Int_t index;
340   for(index = 0; index <fEmcLast-fEmcFirst; index ++)
341     emcExist[index] = 1 ;
342   
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 ;
348   
349   
350   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
351   TIter nextUp(fLinkUpArray) ;
352   
353   AliPHOSLink * linkUp ;
354   
355   AliPHOSCpvRecPoint * nullpointer = 0 ;
356   
357   while ( (linkUp =  static_cast<AliPHOSLink *>(nextUp()) ) ){  
358
359     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
360
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()) ;
366          
367        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
368        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
373     } 
374   }        
375
376   //look through emc recPoints left without CPV
377   if(emcExist){ //if there is emc rec point
378     Int_t iEmcRP ;
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)), 
383                            nullpointer) ;
384        (dynamic_cast<AliPHOSTrackSegment *>(trackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
385        fNTrackSegments++;    
386       } 
387     }
388   }
389   delete [] emcExist ; 
390   delete [] cpvExist ; 
391 }
392
393 //____________________________________________________________________________
394 void  AliPHOSTrackSegmentMakerv1::Exec(Option_t *option)
395 {
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.
400   
401   if(strstr(option,"tim"))
402     gBenchmark->Start("PHOSTSMaker");
403  
404   if(strstr(option,"print")) {
405     Print() ; 
406     return ; 
407   }
408   
409   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
410  
411   const AliPHOSGeometry * geom = gime->PHOSGeometry() ; 
412
413   if (fLastEvent == -1) 
414     fLastEvent = gime->MaxEvent() - 1 ;
415   else 
416     fLastEvent = TMath::Min(fFirstEvent,gime->MaxEvent());
417   Int_t nEvents   = fLastEvent - fFirstEvent + 1;
418
419   Int_t ievent ; 
420   for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
421     gime->Event(ievent,"R") ;
422    //Make some initializations 
423     fNTrackSegments = 0 ;
424     fEmcFirst = 0 ;    
425     fEmcLast  = 0 ;   
426     fCpvFirst = 0 ;   
427     fCpvLast  = 0 ;   
428     
429     gime->TrackSegments()->Clear();
430
431     //    if(!ReadRecPoints(ievent))   continue; //reads RecPoints for event ievent
432     
433     for(fModule = 1; fModule <= geom->GetNModules() ; fModule++ ) {
434       FillOneModule() ; 
435       MakeLinks() ;
436       MakePairs() ;
437     }
438
439     WriteTrackSegments() ;
440
441     if(strstr(option,"deb"))
442       PrintTrackSegments(option);
443     
444     //increment the total number of track segments per run 
445     fTrackSegmentsInRun += gime->TrackSegments()->GetEntriesFast() ; 
446   }
447   
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) ;
453    }
454   if(fWrite) //do not unload in "on flight" mode
455     Unload();
456 }
457
458 //____________________________________________________________________________
459 void AliPHOSTrackSegmentMakerv1::Unload() 
460 {
461   // Unloads the task from the folder
462   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;  
463   gime->PhosLoader()->UnloadRecPoints() ;
464   gime->PhosLoader()->UnloadTracks() ;
465 }
466
467 //____________________________________________________________________________
468 void AliPHOSTrackSegmentMakerv1::Print(const Option_t *)const
469 {
470   //  Print TrackSegmentMaker parameters
471
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) ;
480   }
481   else
482     Info("Print", "AliPHOSTrackSegmentMakerv1 not initialized ") ;
483 }
484
485 //____________________________________________________________________________
486 void AliPHOSTrackSegmentMakerv1::WriteTrackSegments()
487 {
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.
495
496   AliPHOSGetter *gime = AliPHOSGetter::Instance() ; 
497
498   TClonesArray * trackSegments = gime->TrackSegments() ; 
499   trackSegments->Expand(trackSegments->GetEntriesFast()) ;
500
501   if(fWrite){ //We write TreeT
502     TTree * treeT = gime->TreeT();
503     
504     //First TS
505     Int_t bufferSize = 32000 ; 
506     TBranch * tsBranch = treeT->Branch("PHOSTS",&trackSegments,bufferSize);
507     tsBranch->Fill() ;  
508     
509     gime->WriteTracks("OVERWRITE");
510     gime->WriteTrackSegmentMaker("OVERWRITE");
511   }
512 }
513
514
515 //____________________________________________________________________________
516 void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
517 {
518   // option deb - prints # of found TrackSegments
519   // option deb all - prints as well indexed of found RecParticles assigned to the TS
520
521   TClonesArray * trackSegments = AliPHOSGetter::Instance()->TrackSegments() ; 
522
523   Info("PrintTrackSegments", "Results from TrackSegmentMaker:") ; 
524   printf("nevent: %d\n", gAlice->GetEvNumber()) ; 
525   printf("        Found %d TrackSegments\n", trackSegments->GetEntriesFast() ); 
526   
527   if(strstr(option,"all")) {  // printing found TS
528     printf("TrackSegment #  EMC RP#  CPV RP#\n") ; 
529     Int_t index;
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() ) ; 
533     }   
534   }
535 }