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