FairRoot
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PixelFitTracks.cxx
Go to the documentation of this file.
1 /********************************************************************************
2  * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 /*
9  * PixelFitTracks.cxx
10  *
11  * Created on: 25.02.2016
12  * Author: R. Karabowicz
13  */
14 
15 #include "PixelFitTracks.h"
16 
17 #include "FairLogger.h"
18 #include "FairRootManager.h"
19 #include "FairRun.h"
20 #include "FairRuntimeDb.h"
21 #include "PixelDigiPar.h"
22 #include "PixelHit.h"
23 #include "PixelTrack.h"
24 
25 #include <TClonesArray.h>
26 #include <TList.h>
27 #include <TMath.h>
28 
30  : PixelFitTracks("Pixel Track Fitter", 0)
31 {}
32 
34  : PixelFitTracks("Pixel Track Fitter", iVerbose)
35 {}
36 
37 PixelFitTracks::PixelFitTracks(const char* name, Int_t iVerbose)
38  : FairTask(name, iVerbose)
39  , fDigiPar(nullptr)
40  , fHits(nullptr)
41  , fTracks(nullptr)
42  , fFitTracks(nullptr)
43  , fTNofEvents(0)
44  , fNHits(0)
45  , fNTracks(0)
46  , fTNofTracks(0)
47  , fNFitTracks(0)
48  , fTNofFitTracks(0)
49 {
50  Reset();
51 }
52 
54 {
55  Reset();
56  delete fDigiPar;
57  if (fFitTracks) {
58  fFitTracks->Delete();
59  delete fFitTracks;
60  }
61 }
62 
63 void PixelFitTracks::Exec(Option_t* /*opt*/)
64 {
65  Reset();
66 
67  fNHits = fHits->GetEntriesFast();
68  fNTracks = fTracks->GetEntriesFast();
69 
70  LOG(debug) << "PixelFitTracks::Exec() EVENT " << fTNofEvents << " with " << fNTracks << " TRACKS";
71 
72  for (Int_t itrack = 0; itrack < fNTracks; itrack++) {
73  PixelTrack* curTrack = static_cast<PixelTrack*>(fTracks->At(itrack));
74 
75  const Int_t nofHits = curTrack->GetNofHits();
76  // Default initialize the arrays
77  Double_t* hitXPos = new Double_t[nofHits]();
78  Double_t* hitYPos = new Double_t[nofHits]();
79  Double_t* hitZPos = new Double_t[nofHits]();
80 
81  for (Int_t ihit = 0; ihit < nofHits; ihit++) {
82  PixelHit* curHit = static_cast<PixelHit*>(fHits->At(curTrack->GetHitIndex(ihit)));
83 
84  // LOG(info) << " HIT[" << curTrack->GetHitIndex(ihit) << "] = (" << curHit->GetX() << " , " <<
85  // curHit->GetY() << " , " << curHit->GetZ() << ")";
86  hitXPos[ihit] = curHit->GetX();
87  hitYPos[ihit] = curHit->GetY();
88  hitZPos[ihit] = curHit->GetZ();
89  }
90 
91  Double_t valX0 = 0., errX0 = 0., valAX = 0., errAX = 0.;
92  Double_t valY0 = 0., errY0 = 0., valAY = 0., errAY = 0.;
93 
94  LinearRegression(nofHits, hitZPos, hitXPos, valX0, errX0, valAX, errAX);
95  LinearRegression(nofHits, hitZPos, hitYPos, valY0, errY0, valAY, errAY);
96 
97  PixelTrack* fitTrack =
98  new ((*fFitTracks)[fNFitTracks]) PixelTrack(valX0, valAX, valY0, valAY, errX0, errAX, errY0, errAY);
99  for (Int_t ihit = 0; ihit < nofHits; ihit++) {
100  fitTrack->AddHitIndex(curTrack->GetHitIndex(ihit));
101  }
102 
103  fNFitTracks += 1;
104 
105  LOG(debug) << "Track params: "
106  << " AX = " << curTrack->GetAX() << " += " << curTrack->GetAXErr() << " X0 = " << curTrack->GetX0()
107  << " += " << curTrack->GetX0Err() << " "
108  << " AY = " << curTrack->GetAY() << " += " << curTrack->GetAYErr() << " Y0 = " << curTrack->GetY0()
109  << " += " << curTrack->GetY0Err();
110  LOG(debug) << "Fitted params: "
111  << " AX = " << valAX << " += " << errAX << " X0 = " << valX0 << " += " << errX0
112  << " "
113  << " AY = " << valAY << " += " << errAY << " Y0 = " << valY0 << " += " << errY0;
114 
115  delete[] hitXPos;
116  delete[] hitYPos;
117  delete[] hitZPos;
118  }
119 
120  fTNofEvents += 1;
121  fTNofTracks += fNTracks;
122  fTNofFitTracks += fNFitTracks;
123 }
124 
125 Double_t PixelFitTracks::LinearRegression(Int_t nval,
126  Double_t xval[],
127  Double_t yval[],
128  Double_t& valA0,
129  Double_t& errA0,
130  Double_t& valA1,
131  Double_t& errA1)
132 {
133  Double_t valN = static_cast<Double_t>(nval);
134  Double_t sumXY = 0.;
135  Double_t sumX = 0.;
136  Double_t sumY = 0.;
137  Double_t sumXX = 0.;
138 
139  for (Int_t ival = 0; ival < nval; ival++) {
140  sumXY += xval[ival] * yval[ival];
141  sumX += xval[ival];
142  sumY += yval[ival];
143  sumXX += xval[ival] * xval[ival];
144  }
145  valA1 = (valN * sumXY - sumX * sumY) / (valN * sumXX - sumX * sumX);
146  valA0 = (sumY - valA1 * sumX) / valN;
147  Double_t sumEE = 0.;
148  for (Int_t ival = 0; ival < nval; ival++) {
149  sumEE += (yval[ival] - valA0 - valA1 * xval[ival]) * (yval[ival] - valA0 - valA1 * xval[ival]);
150  }
151  Double_t valS = TMath::Sqrt(sumEE / (valN - 2.));
152  errA1 = valS * TMath::Sqrt(valN / (valN * sumXX - sumX * sumX));
153  errA0 = valS * TMath::Sqrt(sumXX / (valN * sumXX - sumX * sumX));
154  // cout << "A0 = " << valA0 << " +- " << errA0 << " / A1 = " << valA1 << " +- " << errA1 << " / S = " << valS <<
155  // endl;
156  return valS;
157 }
158 
159 void PixelFitTracks::SetParContainers()
160 {
161  // Get run and runtime database
162  FairRun* run = FairRun::Instance();
163  if (!run)
164  LOG(fatal) << "No analysis run";
165 
166  FairRuntimeDb* db = run->GetRuntimeDb();
167  if (!db)
168  LOG(fatal) << "No runtime database";
169 
170  // Get GEM digitisation parameter container
171  fDigiPar = static_cast<PixelDigiPar*>(db->getContainer("PixelDigiParameters"));
172 }
173 
174 void PixelFitTracks::GetParList(TList* tempList)
175 {
176  fDigiPar = new PixelDigiPar("PixelDigiParameters");
177  tempList->Add(fDigiPar);
178 
179  return;
180 }
181 
182 void PixelFitTracks::InitMQ(TList* tempList)
183 {
184  LOG(info) << "********************************************** PixelFitTracks::InitMQ()";
185  fDigiPar = (PixelDigiPar*)tempList->FindObject("PixelDigiParameters");
186 
187  fFitTracks = new TClonesArray("PixelTrack", 10000);
188  fFitTracks->SetName("PixelFitTracks");
189  return;
190 }
191 
192 void PixelFitTracks::ExecMQ(TList* inputList, TList* outputList)
193 {
194  // LOG(info) << "********************************************** PixelFitTracks::ExecMQ(" << inputList->GetName() <<
195  // "," << outputList->GetName() << "), Event " << fTNofEvents; LOG(info) <<
196  // "********************************************** PixelFitTracks::ExecMQ(), Event " << fTNofEvents; LOG(info) <<
197  // "f" << FairLogger::flush;
198  fHits = (TClonesArray*)inputList->FindObject("PixelHits");
199  fTracks = (TClonesArray*)inputList->FindObject("PixelTracks");
200  outputList->Add(fFitTracks);
201  Exec("");
202  return;
203 }
204 
205 InitStatus PixelFitTracks::Init()
206 {
207  // Get input array
209  if (!ioman)
210  LOG(fatal) << "No FairRootManager";
211 
212  fHits = static_cast<TClonesArray*>(ioman->GetObject("PixelHits"));
213  if (!fHits)
214  LOG(warn) << "PixelFitTracks::Init() No input PixelHit array!";
215  fTracks = static_cast<TClonesArray*>(ioman->GetObject("PixelTracks"));
216  if (!fTracks)
217  LOG(warn) << "PixelFitTracks::Init() No input PixelTrack array!";
218 
219  // Register output array PixelHit
220  fFitTracks = new TClonesArray("PixelTrack", 10000);
221  ioman->Register("PixelFitTracks", "Pixel", fFitTracks, kTRUE);
222 
223  return kSUCCESS;
224 }
225 
226 InitStatus PixelFitTracks::ReInit() { return kSUCCESS; }
227 
228 void PixelFitTracks::Reset()
229 {
230  fNFitTracks = fNTracks = fNHits = 0;
231  if (fFitTracks)
232  fFitTracks->Clear();
233 }
234 
235 void PixelFitTracks::Finish()
236 {
237  if (fFitTracks)
238  fFitTracks->Delete();
239 
240  LOG(info) << "-------------------- " << fName.Data() << " : Summary ------------------------";
241  LOG(info) << " Events: " << fTNofEvents;
242  LOG(info) << " Tracks: " << fTNofTracks << " ("
243  << static_cast<Double_t>(fTNofTracks) / (static_cast<Double_t>(fTNofEvents)) << " per event)";
244  LOG(info) << " Fitted Tracks: " << fTNofFitTracks << " ("
245  << static_cast<Double_t>(fTNofFitTracks) / (static_cast<Double_t>(fTNofEvents)) << " per event)";
246  LOG(info) << "---------------------------------------------------------------------";
247 }
248 
list of container factories
Definition: FairRuntimeDb.h:24
Double_t GetZ() const
Definition: FairHit.h:50
void AddHitIndex(Int_t hitIndex)
Definition: PixelTrack.h:64
InitStatus
Definition: FairTask.h:33
Double_t GetX0Err()
Definition: PixelTrack.h:59
static FairRun * Instance()
Definition: FairRun.cxx:31
static FairRootManager * Instance()
ClassImp(FairEventBuilder)
virtual ~PixelFitTracks()
Double_t GetX() const
Definition: FairHit.h:48
Digitization Parameter Class for Pixel detector.
Definition: PixelDigiPar.h:24
TObject * GetObject(const char *BrName)
FairParSet * getContainer(const Text_t *)
Int_t GetHitIndex(Int_t ihit)
Definition: PixelTrack.h:66
Int_t GetNofHits()
Definition: PixelTrack.h:65
virtual void Exec(Option_t *opt)
Double_t GetAX()
Definition: PixelTrack.h:56
virtual void ExecMQ(TList *inputList, TList *outputList)
FairRuntimeDb * GetRuntimeDb(void)
Definition: FairRun.h:80
Double_t GetY() const
Definition: FairHit.h:49
Double_t GetY0()
Definition: PixelTrack.h:57
void Register(const char *name, const char *Foldername, TNamed *obj, Bool_t toFile)
Double_t GetAY()
Definition: PixelTrack.h:58
virtual void InitMQ(TList *tempList)
Double_t GetX0()
Definition: PixelTrack.h:55
Double_t GetAXErr()
Definition: PixelTrack.h:60
Double_t GetY0Err()
Definition: PixelTrack.h:61
virtual void GetParList(TList *tempList)
Double_t GetAYErr()
Definition: PixelTrack.h:62