Fix bugs in PID assignment
[u/mrichter/AliRoot.git] / TRD / AliTRDReconstructor.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
16 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 // class for TRD reconstruction                                              //
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include <TFile.h>
25
26 #include "AliTRDReconstructor.h"
27 #include "AliRunLoader.h"
28 #include "AliTRDclusterizerV1.h"
29 #include "AliTRDtracker.h"
30 #include "AliTRDpidESD.h"
31 #include "AliRawReader.h"
32 #include "AliLog.h"
33 #include "AliTRDtrigger.h"
34 #include "AliTRDtrigParam.h"
35 #include "AliTRDgtuTrack.h"
36 #include "AliRun.h"
37 #include "AliESDTrdTrack.h"
38 #include "AliESD.h"
39
40 ClassImp(AliTRDReconstructor)
41
42 Bool_t  AliTRDReconstructor::fgkSeedingOn = kFALSE;
43 Int_t   AliTRDReconstructor::fgStreamLevel     = 0;      // stream (debug) level
44
45 //_____________________________________________________________________________
46 void AliTRDReconstructor::Reconstruct(AliRunLoader* runLoader) const
47 {
48 // reconstruct clusters
49
50   AliLoader *loader=runLoader->GetLoader("TRDLoader");
51   loader->LoadRecPoints("recreate");
52
53   runLoader->CdGAFile();
54   Int_t nEvents = runLoader->GetNumberOfEvents();
55
56   for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
57     AliTRDclusterizerV1 clusterer("clusterer", "TRD clusterizer");
58     clusterer.Open(runLoader->GetFileName(), iEvent);
59     clusterer.ReadDigits();
60     clusterer.MakeClusters();
61     clusterer.WriteClusters(-1);
62   }
63
64   loader->UnloadRecPoints();
65
66   // Trigger (tracklets, LTU)
67
68   loader->LoadTracks("RECREATE");
69
70   Info("Reconstruct","Trigger tracklets will be produced");
71
72   AliTRDtrigger trdTrigger("Trigger","Trigger class"); 
73
74   AliTRDtrigParam *trigp = new AliTRDtrigParam("TRDtrigParam","TRD Trigger parameters");
75
76   if (runLoader->GetAliRun() == 0x0) runLoader->LoadgAlice();
77   gAlice = runLoader->GetAliRun();
78   Double_t x[3] = { 0.0, 0.0, 0.0 };
79   Double_t b[3];
80   gAlice->Field(x,b);  // b[] is in kilo Gauss
81   Float_t field = b[2] * 0.1; // Tesla
82   Info("Reconstruct","Trigger set for magnetic field = %f Tesla \n",field);
83
84   trigp->SetField(field);
85   trigp->Init();
86   trdTrigger.SetParameter(trigp);
87
88   for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
89     trdTrigger.Open(runLoader->GetFileName(), iEvent);
90     trdTrigger.ReadDigits();
91     trdTrigger.MakeTracklets();
92     trdTrigger.WriteTracklets(-1);
93   }
94
95   loader->UnloadTracks();
96
97 }
98
99 //_____________________________________________________________________________
100 void AliTRDReconstructor::Reconstruct(AliRunLoader* runLoader,
101                                       AliRawReader* rawReader) const
102 {
103 // reconstruct clusters
104
105   AliInfo("Reconstruct TRD clusters from RAW data");
106
107   AliLoader *loader=runLoader->GetLoader("TRDLoader");
108   loader->LoadRecPoints("recreate");
109
110   runLoader->CdGAFile();
111   Int_t nEvents = runLoader->GetNumberOfEvents();
112
113   for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
114     if (!rawReader->NextEvent()) break;
115     AliTRDclusterizerV1 clusterer("clusterer", "TRD clusterizer");
116     clusterer.Open(runLoader->GetFileName(), iEvent);
117     clusterer.ReadDigits(rawReader);
118     clusterer.MakeClusters();
119     clusterer.WriteClusters(-1);
120   }
121
122   loader->UnloadRecPoints();
123
124   // Trigger (tracklets, LTU)
125
126   loader->LoadTracks();
127   if (loader->TreeT()) {
128     Info("Reconstruct","Tracklets already exist");
129     return;
130   }
131   Info("Reconstruct","Trigger tracklets will be produced");
132
133   AliTRDtrigger trdTrigger("Trigger","Trigger class"); 
134
135   AliTRDtrigParam *trigp = new AliTRDtrigParam("TRDtrigParam","TRD Trigger parameters");
136
137   if (runLoader->GetAliRun() == 0x0) runLoader->LoadgAlice();
138   gAlice = runLoader->GetAliRun();
139   Double_t x[3] = { 0.0, 0.0, 0.0 };
140   Double_t b[3];
141   gAlice->Field(x,b);  // b[] is in kilo Gauss
142   Float_t field = b[2] * 0.1; // Tesla
143   Info("Reconstruct","Trigger set for magnetic field = %f Tesla \n",field);
144
145   trigp->SetField(field);
146   trigp->Init();
147   trdTrigger.SetParameter(trigp);
148
149   rawReader->RewindEvents();
150
151   for (Int_t iEvent = 0; iEvent < nEvents; iEvent++) {
152     if (!rawReader->NextEvent()) break;
153     trdTrigger.Open(runLoader->GetFileName(), iEvent);
154     trdTrigger.ReadDigits(rawReader);
155     trdTrigger.MakeTracklets();
156     trdTrigger.WriteTracklets(-1);
157   }
158
159   loader->UnloadTracks();
160
161 }
162
163 //_____________________________________________________________________________
164 AliTracker* AliTRDReconstructor::CreateTracker(AliRunLoader* runLoader) const
165 {
166 // create a TRD tracker
167
168   runLoader->CdGAFile();
169   return new AliTRDtracker(gFile);
170 }
171
172 //_____________________________________________________________________________
173 void AliTRDReconstructor::FillESD(AliRunLoader* runLoader, 
174                                   AliESD* esd) const
175 {
176 // make PID
177
178   //Double_t parTRD[] = {
179   //  280., // Min. Ionizing Particle signal.  Check it !!!
180   //  0.23, // relative resolution             Check it !!!
181   //  10.   // PID range (in sigmas)
182   //};
183   AliTRDpidESD trdPID;
184   trdPID.MakePID(esd);
185
186   // Trigger (tracks, GTU)
187
188   AliTRDtrigger trdTrigger("Trigger","Trigger class"); 
189
190   AliTRDtrigParam *trigp = new AliTRDtrigParam("TRDtrigParam","TRD Trigger parameters");
191
192   if (runLoader->GetAliRun() == 0x0) runLoader->LoadgAlice();
193   gAlice = runLoader->GetAliRun();
194   Double_t x[3] = { 0.0, 0.0, 0.0 };
195   Double_t b[3];
196   gAlice->Field(x,b);  // b[] is in kilo Gauss
197   Float_t field = b[2] * 0.1; // Tesla
198   Info("FillESD","Trigger set for magnetic field = %f Tesla \n",field);
199
200   trigp->SetField(field);
201   trigp->Init();
202
203   trdTrigger.SetParameter(trigp);
204   trdTrigger.SetRunLoader(runLoader);
205   trdTrigger.Init();
206
207   Int_t iEvent = runLoader->GetEventNumber(); 
208   runLoader->GetEvent(iEvent);
209   trdTrigger.ReadTracklets(runLoader);
210
211   AliESDTrdTrack *TrdTrack = new AliESDTrdTrack();
212   AliTRDgtuTrack *GtuTrack;
213
214   Int_t nTracks = trdTrigger.GetNumberOfTracks();
215   for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
216
217     GtuTrack = trdTrigger.GetTrack(iTrack);
218
219     TrdTrack->SetYproj(GtuTrack->GetYproj());
220     TrdTrack->SetZproj(GtuTrack->GetZproj());
221     TrdTrack->SetSlope(GtuTrack->GetSlope());
222     TrdTrack->SetDetector(GtuTrack->GetDetector());
223     TrdTrack->SetTracklets(GtuTrack->GetTracklets());
224     TrdTrack->SetPlanes(GtuTrack->GetPlanes());
225     TrdTrack->SetClusters(GtuTrack->GetClusters());
226     TrdTrack->SetPt(GtuTrack->GetPt());
227     TrdTrack->SetPhi(GtuTrack->GetPhi());
228     TrdTrack->SetEta(GtuTrack->GetEta());
229     TrdTrack->SetLabel(GtuTrack->GetLabel());
230     TrdTrack->SetPID(GtuTrack->GetPID());
231     TrdTrack->SetIsElectron(GtuTrack->IsElectron());
232
233     esd->AddTrdTrack(TrdTrack);
234
235   }
236
237 }
238
239
240