0951c096063bc73b33ffdbe2e05694f885b6e5d7
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTracker.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.5  2007/08/03 13:52:16  kharlov
21  * Working skeleton of matching the ESD tracks and ESD clusters (Iouri Belikov)
22  *
23  */
24
25 #include <TClonesArray.h>
26 #include <TMath.h>
27
28 #include <AliLog.h>
29 #include "AliPHOSTracker.h"
30 #include "AliPHOSEmcRecPoint.h"
31 #include "AliESDEvent.h"
32 #include "AliESDtrack.h"
33
34 //-------------------------------------------------------------------------
35 //                          PHOS tracker.
36 // Matches ESD tracks with the PHOS and makes the PID.  
37 //
38 //-------------------------------------------------------------------------
39
40 ClassImp(AliPHOSTracker)
41
42 Bool_t AliPHOSTracker::fgDebug = kFALSE ;  
43
44
45 // ***** Some geometrical constants (used in PropagateBack) 
46
47 const Double_t kR=460.+ 9;  // Radial coord. of the centre of EMC module (cm)
48
49 const Double_t kAlpha=20.*TMath::Pi()/180.;     // Segmentation angle (rad)
50 const Double_t kYmax=kR*TMath::Tan(0.5*kAlpha); // Maximal possible y-coord.(cm)
51 const Double_t kZmax=65.; // Approximately: the maximal possible z-coord.(cm)
52
53
54
55 AliPHOSTracker::AliPHOSTracker(): AliTracker(), fRunLoader(0) {
56   //--------------------------------------------------------------------
57   // The default constructor
58   //--------------------------------------------------------------------
59   for (Int_t i=0; i<5; i++) 
60       fModules[i]=new TClonesArray("AliPHOSEmcRecPoint",777);
61 }
62
63 AliPHOSTracker::~AliPHOSTracker() {
64   //--------------------------------------------------------------------
65   // The destructor
66   //--------------------------------------------------------------------
67   for (Int_t i=0; i<5; i++) {
68       (fModules[i])->Delete();
69       delete fModules[i];
70   }
71 }
72
73 Int_t AliPHOSTracker::LoadClusters(TTree *cTree) {
74   //--------------------------------------------------------------------
75   // This function loads the PHOS clusters
76   //--------------------------------------------------------------------
77   TObjArray *arr=NULL;
78   TBranch *branch=cTree->GetBranch("PHOSEmcRP");
79   if (branch==0) {
80     AliError("No branch with the EMC clusters found !");
81     return 1;
82   }
83   branch->SetAddress(&arr);
84
85   Int_t nclusters=0;
86   Int_t nentr=(Int_t)branch->GetEntries();
87   for (Int_t i=0; i<nentr; i++) {
88     if (!branch->GetEvent(i)) continue;
89     Int_t ncl=arr->GetEntriesFast();
90     while (ncl--) {
91       AliPHOSEmcRecPoint *cl=(AliPHOSEmcRecPoint*)arr->UncheckedAt(ncl);
92
93       Int_t m=cl->GetPHOSMod();
94       if ((m<1)||(m>5)) {
95          AliError("Wrong module index !");
96          return 1;
97       }
98
99       // Here is how the alignment is treated
100       if (!cl->Misalign()) AliWarning("Can't misalign this cluster !");
101
102       cl->SetBit(14,kFALSE); // The clusters are not yet attached to any track
103
104       TClonesArray &module=*fModules[m-1];
105       Int_t idx=module.GetEntriesFast();
106       new (module[idx]) AliPHOSEmcRecPoint(*cl); 
107
108       nclusters++;
109
110     }
111   }  
112
113   Info("LoadClusters","Number of loaded clusters: %d",nclusters);
114
115   return 0;
116 }
117
118
119 Int_t AliPHOSTracker::PropagateBack(AliESDEvent *esd) {
120   //--------------------------------------------------------------------
121   // Called by AliReconstruction 
122   // Performs the track matching with the PHOS modules
123   // Makes the PID
124   //--------------------------------------------------------------------
125
126   // The following old function is a bad function: it uses RunLoader ;(
127   PropagateBackOld(esd);   
128
129
130   Int_t nt=esd->GetNumberOfTracks();
131
132   // *** Select and sort the ESD track in accordance with their quality
133   Double_t *quality=new Double_t[nt];
134   Int_t *index=new Int_t[nt];  
135   for (Int_t i=0; i<nt; i++) {
136      AliESDtrack *esdTrack=esd->GetTrack(i);
137      quality[i] = esdTrack->GetSigmaY2() + esdTrack->GetSigmaZ2();
138   }
139   TMath::Sort(nt,quality,index,kFALSE);
140
141
142   // *** Start the matching
143   Double_t bz=GetBz(); 
144   Int_t matched=0;
145   for (Int_t i=0; i<nt; i++) {
146      AliESDtrack *esdTrack=esd->GetTrack(index[i]);
147
148      // Skip the tracks having "wrong" status (has to be checked/tuned)
149      ULong_t status = esdTrack->GetStatus();
150      if ((status & AliESDtrack::kTRDout)   == 0) continue;
151      if ((status & AliESDtrack::kTRDrefit) == 1) continue;
152
153      AliExternalTrackParam t(*esdTrack);
154
155      Int_t isec=Int_t(t.GetAlpha()/kAlpha);
156      Int_t imod=-isec-2; // PHOS module
157
158      Double_t y;                       // Some tracks do not reach the PHOS
159      if (!t.GetYAt(kR,bz,y)) continue; //    because of the bending
160
161      Double_t z; t.GetZAt(kR,bz,z);   
162      if (TMath::Abs(z) > kZmax) continue; // Some tracks miss the PHOS in Z
163  
164      Bool_t ok=kTRUE;
165      while (TMath::Abs(y) > kYmax) {   // Find the matching module
166         Double_t alp=t.GetAlpha();
167         if (y > kYmax) {
168           if (!t.Rotate(alp+kAlpha)) {ok=kFALSE; break;}
169           imod--;
170         } else if (y < -kYmax) {
171           if (!t.Rotate(alp-kAlpha)) {ok=kFALSE; break;}
172           imod++;
173         }
174         if (!t.GetYAt(kR,bz,y)) {ok=kFALSE; break;}
175      }
176      if (!ok) continue; // Track rotation failed
177  
178   
179      if ((imod<0)||(imod>4)) continue; // Some tracks miss the PHOS in azimuth
180
181      //t.CorrectForMaterial(...); // Correct for the TOF material, if needed
182      t.PropagateTo(kR,bz);        // Propagate to the matching module
183
184
185     // *** Search for the "best" cluster (can be improved)
186      TClonesArray &cArray=*fModules[imod];
187      Int_t ncl=cArray.GetEntriesFast();
188      AliPHOSEmcRecPoint *bestCluster=0;            // The "best" cluster
189      Double_t maxd2=400; // (cm^2)
190      for (Int_t i=0; i<ncl; i++) {
191        AliPHOSEmcRecPoint *c=(AliPHOSEmcRecPoint *)cArray.UncheckedAt(i);
192
193        if (c->TestBit(14)) continue; // This clusters is "used"
194
195        Double_t dy = t.GetY() - c->GetY(), dz = t.GetZ() - c->GetZ();
196        Double_t d2 = dy*dy + dz*dz;
197        if (d2 < maxd2) {
198           maxd2=d2;
199           bestCluster=c;
200        }
201      }
202
203      if (!bestCluster) continue;   // No reasonable matching found 
204
205      bestCluster->SetBit(14,kTRUE); // This clusters is now attached to a track
206
207      matched++;
208
209      // *** Now, do the PID with the "bestCluster"
210      // and add the corresponding info to the ESD track pointed by "esdTrack"  
211
212      /*
213      printf("%e %e %e %e\n",t.GetSign(), t.GetX() - bestCluster->GetX(),
214                                          t.GetY() - bestCluster->GetY(),
215                                          t.GetZ() - bestCluster->GetZ());
216      */
217   }
218     
219   Info("PropagateBack","Number of matched tracks: %d",matched);
220
221   delete[] quality;
222   delete[] index;
223
224   return 0;
225 }
226
227 AliCluster *AliPHOSTracker::GetCluster(Int_t index) const {
228   //--------------------------------------------------------------------
229   // Returns the pointer to a given cluster
230   //--------------------------------------------------------------------
231   Int_t m=(index & 0xf0000000) >> 28;  // Module number
232   Int_t i=(index & 0x0fffffff) >> 00;  // Index within the module
233   
234   return (AliCluster*)(fModules[m])->UncheckedAt(i);
235 }
236
237 void AliPHOSTracker::UnloadClusters() {
238   //--------------------------------------------------------------------
239   // This function unloads the PHOS clusters
240   //--------------------------------------------------------------------
241   for (Int_t i=0; i<5; i++) (fModules[i])->Delete();
242 }
243
244
245
246 // **** The following are bad functions:  they use RunLoader ;(
247 // **** To be rewritten.
248
249 #include "AliPHOSTrackSegmentMakerv1.h"
250 #include "AliPHOSTrackSegmentMakerv2.h"
251 #include "AliPHOSPIDv1.h"
252 #include "AliRunLoader.h"
253
254 Int_t AliPHOSTracker::PropagateBackOld(AliESDEvent *esd) {
255   // Bad function: it uses RunLoader ;(  
256   // Creates the tracksegments and Recparticles
257   // Makes the PID
258   
259   Int_t eventNumber = fRunLoader->GetEventNumber() ;
260
261   TString headerFile(fRunLoader->GetFileName()) ; 
262   TString branchName(fRunLoader->GetEventFolder()->GetName()) ;  
263  
264   AliPHOSTrackSegmentMakerv1 tsm(headerFile, branchName);
265 //  AliPHOSTrackSegmentMakerv2 tsm(headerFile, branchName);
266   tsm.SetESD(esd) ; 
267   AliPHOSPIDv1 pid(headerFile, branchName);
268   pid.SetESD(esd) ; 
269
270   //PH  SetDebug() ;
271
272   // do current event; the loop over events is done by AliReconstruction::Run()
273   tsm.SetEventRange(eventNumber, eventNumber) ; 
274   pid.SetEventRange(eventNumber, eventNumber) ; 
275   if ( Debug() ) {
276    tsm.ExecuteTask("deb all") ;
277    pid.ExecuteTask("deb all") ;
278   }
279   else {
280     tsm.ExecuteTask("") ;
281     pid.ExecuteTask("") ;
282   }
283   
284   return 0;
285 }
286
287 AliPHOSTracker::AliPHOSTracker(AliRunLoader *l): AliTracker(), fRunLoader(l) {
288   //--------------------------------------------------------------------
289   // Bad constructor:  uses RunLoader ;(
290   //--------------------------------------------------------------------
291   for (Int_t i=0; i<5; i++) 
292       fModules[i]=new TClonesArray("AliPHOSEmcRecPoint",777);
293 }
294