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