]>
Commit | Line | Data |
---|---|---|
59fbd3a7 | 1 | #define PointMaker_cxx |
2 | // The class definition in esdTree.h has been generated automatically | |
3 | // by the ROOT utility TTree::MakeSelector(). This class is derived | |
4 | // from the ROOT class TSelector. For more information on the TSelector | |
5 | // framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual. | |
6 | ||
7 | // The following methods are defined in this file: | |
8 | // Begin(): called everytime a loop on the tree starts, | |
9 | // a convenient place to create your histograms. | |
10 | // SlaveBegin(): called after Begin(), when on PROOF called only on the | |
11 | // slave servers. | |
12 | // Process(): called for each event, in this function you decide what | |
13 | // to read and fill your histograms. | |
14 | // SlaveTerminate: called at the end of the loop on the tree, when on PROOF | |
15 | // called only on the slave servers. | |
16 | // Terminate(): called at the end of the loop on the tree, | |
17 | // a convenient place to draw/fit your histograms. | |
18 | // | |
19 | ||
20 | // To do: rewrite as AliAnalysisTaskSE.(like initiated by Alex at some point) | |
21 | ||
22 | #include <TROOT.h> | |
23 | #include <TStyle.h> | |
24 | #include <TLinearFitter.h> | |
25 | #include "AliAlignObjParams.h" | |
26 | #include "AliTrackPointArray.h" | |
27 | #include "AliLog.h" | |
28 | #include "PointMaker.h" | |
29 | //============================================================================= | |
30 | PointMaker::PointMaker(char *outfil) : | |
31 | TSelector(), | |
32 | fChain(0), | |
33 | fESD(0), | |
34 | // fESDfriend(0), | |
35 | fFile(0), | |
36 | fTree(0), | |
37 | fArray(0), | |
38 | fNevents(0), | |
39 | fNtracks(0), | |
40 | fNAcceptedTracks(0), | |
41 | fOutfil(outfil) | |
42 | { | |
43 | // Constructor. Initialization of pointers | |
44 | } | |
45 | //============================================================================= | |
46 | PointMaker::~PointMaker() { | |
47 | // Remove all pointers | |
48 | ||
49 | // delete fESD; | |
50 | ||
51 | // histograms are in the output list and deleted when the output | |
52 | // list is deleted by the TSelector dtor | |
53 | delete fFile; | |
54 | } | |
55 | //============================================================================= | |
56 | void PointMaker::Begin(TTree *) | |
57 | { | |
58 | // The Begin() function is called at the start of the query. | |
59 | // When running with PROOF Begin() is only called on the client. | |
60 | // The tree argument is deprecated (on PROOF 0 is passed). | |
61 | ||
62 | TString option = GetOption(); | |
63 | ||
64 | // create output file and tree - TEMPORARY -until root trees are merged in the memory | |
65 | fFile = TFile::Open(fOutfil.Data(), "RECREATE",0); | |
66 | fTree = new TTree("spTree", "Tree with track space point arrays"); | |
67 | fTree->Branch("SP","AliTrackPointArray", &fArray); | |
68 | printf("Begin called\n"); | |
69 | } | |
70 | //============================================================================= | |
71 | void PointMaker::SlaveBegin(TTree * tree) | |
72 | { | |
73 | // The SlaveBegin() function is called after the Begin() function. | |
74 | // When running with PROOF SlaveBegin() is called on each slave server. | |
75 | // The tree argument is deprecated (on PROOF 0 is passed). | |
76 | ||
77 | Init(tree); | |
78 | ||
79 | AliAlignObjParams alobj; // initialize align obj. | |
80 | TString option = GetOption(); | |
81 | ||
82 | // histograms to monitor cut efficiency and module population | |
83 | fCuttra = new TH1D("cuttra","cuttra",20,0.5,20.5); | |
84 | fCutpoi = new TH1D("cutpoi","cutpoi",20,0.5,20.5); | |
85 | fModpop = new TH2D("modpop","modpop",90,-0.5,89.5,30,-0.5,29.5); | |
86 | fModpop->SetXTitle("module nr"); | |
87 | fModpop->SetYTitle("layer nr"); | |
88 | printf("SlaveBegin called\n"); | |
89 | } | |
90 | //============================================================================= | |
91 | void PointMaker::Init(TTree *tree) | |
92 | { | |
93 | // The Init() function is called when the selector needs to initialize | |
94 | // a new tree or chain. Typically here the branch addresses of the tree | |
95 | // will be set. It is normaly not necessary to make changes to the | |
96 | // generated code, but the routine can be extended by the user if needed. | |
97 | // Init() will be called many times when running with PROOF. | |
98 | ||
99 | // Set branch addresses | |
100 | if (tree == 0) return; | |
101 | fChain = tree; | |
102 | ||
103 | /* | |
104 | fChain->SetBranchAddress("ESD",&fESD); | |
105 | fChain->SetBranchStatus("ESDfriend*",1); | |
106 | fChain->SetBranchAddress("ESDfriend.",&fESDfriend); | |
107 | */ | |
108 | ||
109 | fESD = new AliESDEvent(); | |
110 | fChain->SetBranchStatus("ESDfriend*",1); | |
111 | fESD->ReadFromTree(fChain); | |
112 | fESDfriend = (AliESDfriend*)fESD->FindListObject("AliESDfriend"); | |
113 | if(!fESDfriend) fChain->SetBranchAddress("ESDfriend.",&fESDfriend); | |
114 | printf("Init called\n"); | |
115 | } | |
116 | //============================================================================= | |
117 | Bool_t PointMaker::Notify() | |
118 | { | |
119 | // The Notify() function is called when a new file is opened. This | |
120 | // can be either for a new TTree in a TChain or when when a new TTree | |
121 | // is started when using PROOF. Typically here the branch pointers | |
122 | // will be retrieved. It is normaly not necessary to make changes | |
123 | // to the generated code, but the routine can be extended by the | |
124 | // user if needed. | |
125 | ||
126 | printf("Notify called\n"); | |
127 | return kTRUE; | |
128 | } | |
129 | //============================================================================= | |
130 | Bool_t PointMaker::IsIdenticalWithOneOf(AliTrackPoint *p, AliTrackPointArray *parray, int nmax) { | |
131 | ||
132 | // Is the point p identical with one of the points on the list parray? | |
133 | // This is a fix for aliroot 4-16-Rev-01 (and before) writing some | |
134 | // spurious uninitialized points. | |
135 | ||
136 | for (int i=0; i<parray->GetNPoints() && i<nmax; i++) { | |
137 | AliTrackPoint pa; | |
138 | parray->GetPoint(pa,i); | |
139 | //printf("comparing %7.3f with %7.3f\n",p->GetY(),pa.GetY()); | |
140 | if (p->GetResidual(pa,0)<1e-8) return kTRUE; | |
141 | //printf("different\n"); | |
142 | } | |
143 | return kFALSE; | |
144 | } | |
145 | //============================================================================= | |
146 | Bool_t PointMaker::Process(Long64_t entry) | |
147 | { | |
148 | // The Process() function is called for each entry in the tree (or possibly | |
149 | // keyed object in the case of PROOF) to be processed. The entry argument | |
150 | // specifies which entry in the currently loaded tree is to be processed. | |
151 | // It can be passed to either TTree::GetEntry() or TBranch::GetEntry() | |
152 | // to read either all or the required parts of the data. When processing | |
153 | // keyed objects with PROOF, the object is already loaded and is available | |
154 | // via the fObject pointer. | |
155 | // | |
156 | // This function should contain the "body" of the analysis. It can contain | |
157 | // simple or elaborate selection criteria, run algorithms on the data | |
158 | // of the event and typically fill histograms. | |
159 | ||
160 | // WARNING when a selector is used with a TChain, you must use | |
161 | // the pointer to the current TTree to call GetEntry(entry). | |
162 | // The entry is always the local entry number in the current tree. | |
163 | // Assuming that fChain is the pointer to the TChain being processed, | |
164 | // use fChain->GetTree()->GetEntry(entry). | |
165 | ||
166 | fChain->GetTree()->GetEntry(entry); | |
167 | fESD->SetESDfriend(fESDfriend); | |
168 | ||
169 | TLinearFitter fitter(2, "pol1"); | |
170 | TLinearFitter fitterz(2, "pol1"); | |
171 | ||
172 | // track cuts | |
173 | ||
174 | int tpc = 0; // require tpc | |
175 | int ptu = 0; // require certain pt's (magnetic field and tpc presumably on) | |
176 | ||
177 | const Float_t kMaxDelta = 1; | |
178 | const Float_t kMinNcl = 60; | |
179 | const Float_t kMinPtLow = 0.2; | |
180 | const Float_t kMinNclLow = 100; | |
181 | const Float_t kMinPt0 = 2; | |
182 | UInt_t status = AliESDtrack::kTRDrefit; | |
183 | if (tpc) status |= AliESDtrack::kTPCrefit; | |
184 | ||
185 | const Float_t kMinRadius2 = 2*2; | |
186 | const Float_t kMaxRadius2 = 400*400; | |
187 | const Float_t kDeadSpace = 4; | |
188 | const Float_t kTan = TMath::Tan(10*TMath::DegToRad()); | |
189 | Int_t ntracks = fESD->GetNumberOfTracks(); | |
190 | const AliTrackPointArray *array = 0; | |
191 | AliTrackPointArray *tmpArray = 0; | |
192 | // trdarray contains all trd points in this event, for duplication detection | |
193 | AliTrackPointArray *trdarray = new AliTrackPointArray(1000); | |
194 | int ntrdarray = 0; | |
195 | ||
196 | int pr = (!gROOT->IsBatch()); | |
197 | for (Int_t itrack=0; itrack < ntracks; itrack++) { | |
198 | if (pr) printf("\revent %d track %d",fNevents,itrack); | |
199 | fCuttra->Fill(1); | |
200 | AliESDtrack * track = fESD->GetTrack(itrack); | |
201 | fNtracks++; | |
202 | if (!track) continue; | |
203 | fCuttra->Fill(2); | |
204 | //if ((track->GetStatus() & status) == 0) continue; | |
205 | fCuttra->Fill(3); | |
206 | if (track->GetKinkIndex(0)!=0) continue; | |
207 | fCuttra->Fill(4); | |
208 | if (tpc) if (track->GetTPCNcls()<kMinNcl) continue; | |
209 | fCuttra->Fill(5); | |
210 | if (ptu) if (track->GetP() < kMinPtLow) continue; | |
211 | fCuttra->Fill(6); | |
212 | if (ptu) if (track->GetP() < kMinPt0 && track->GetTPCNcls()<kMinNclLow) continue; | |
213 | fCuttra->Fill(7); | |
214 | // | |
215 | // select points | |
216 | // | |
217 | array = track->GetTrackPointArray(); | |
218 | if (!array) continue; | |
219 | Int_t npoints = array->GetNPoints(); | |
220 | if (tmpArray) delete tmpArray; | |
221 | tmpArray = new AliTrackPointArray(npoints); | |
222 | Int_t current = 0; | |
223 | int ntpc = 0; // number of good TPC points | |
224 | int ntrd = 0; // number of good TRD points | |
225 | for (Int_t ipoint=0; ipoint<npoints; ipoint++){ | |
226 | fCutpoi->Fill(1); | |
227 | AliTrackPoint p; | |
228 | array->GetPoint(p, ipoint); | |
229 | ||
230 | // filter buggy points | |
231 | ||
232 | UShort_t volid = array->GetVolumeID()[ipoint]; | |
233 | Int_t iModule; | |
234 | AliGeomManager::ELayerID layer = AliGeomManager::VolUIDToLayer(volid,iModule); | |
235 | if ((layer < AliGeomManager::kFirstLayer) || (layer >= AliGeomManager::kLastLayer)) continue; | |
236 | fCutpoi->Fill(2); | |
237 | if ((iModule >= AliGeomManager::LayerSize(layer)) || (iModule < 0)) continue; | |
238 | fCutpoi->Fill(3); | |
239 | Float_t r2 = p.GetX()*p.GetX()+p.GetY()*p.GetY(); | |
240 | if ( r2<kMinRadius2 || r2 > kMaxRadius2 ) continue; | |
241 | fCutpoi->Fill(4); | |
242 | ||
243 | // TPC clusters handling | |
244 | ||
245 | if (layer>=AliGeomManager::kTPC1 && layer<=AliGeomManager::kTPC2){ | |
246 | // invalid covariance | |
247 | if (p.GetCov()[0]<0 || p.GetCov()[3]<0 || p.GetCov()[5]<0) continue; | |
248 | fCutpoi->Fill(7); | |
249 | ||
250 | // remove edge clusters | |
251 | ||
252 | AliTrackPoint& plocal = p.MasterToLocal(); | |
253 | Double_t ylocal = plocal.GetY(); | |
254 | Double_t zlocal = plocal.GetZ(); | |
255 | Double_t xlocal = plocal.GetX(); | |
256 | Float_t edgey = TMath::Abs(plocal.GetX()*kTan); | |
257 | Int_t nclose=0; | |
258 | fitter.ClearPoints(); | |
259 | fitterz.ClearPoints(); | |
260 | for (Int_t jpoint=ipoint-7; jpoint<=ipoint+7; jpoint++){ | |
261 | if (jpoint<0 || jpoint>=npoints) continue; | |
262 | if (ipoint==jpoint) continue; | |
263 | UShort_t volidL = array->GetVolumeID()[jpoint]; | |
264 | if (volidL!=volid) continue; | |
265 | AliTrackPoint pc; | |
266 | array->GetPoint(pc, jpoint); | |
267 | AliTrackPoint &pcl= pc.MasterToLocal(); | |
268 | Double_t dx = pcl.GetX()-xlocal; | |
269 | fitter.AddPoint(&dx,pcl.GetY(),1); | |
270 | fitterz.AddPoint(&dx,pcl.GetZ(),1); | |
271 | nclose++; | |
272 | } | |
273 | if (nclose<6) continue; | |
274 | fCutpoi->Fill(8); | |
275 | fitter.Eval(); | |
276 | fitterz.Eval(); | |
277 | Double_t fity =fitter.GetParameter(0); | |
278 | Double_t fitz =fitterz.GetParameter(0); | |
279 | if (TMath::Abs(ylocal-fity)>kMaxDelta) continue; //outlier | |
280 | fCutpoi->Fill(9); | |
281 | if (TMath::Abs(zlocal-fitz)>kMaxDelta) continue; //outlier | |
282 | fCutpoi->Fill(10); | |
283 | if (TMath::Abs(fity)>edgey-kDeadSpace) continue; // remove edge clusters | |
284 | fCutpoi->Fill(11); | |
285 | ntpc++; | |
286 | } | |
287 | fCutpoi->Fill(12); | |
288 | ||
289 | // TRD track points handling | |
290 | ||
291 | if (layer>=AliGeomManager::kTRD1 && layer<=AliGeomManager::kTRD6){ | |
292 | fCutpoi->Fill(14); | |
293 | if (IsIdenticalWithOneOf(&p,trdarray,ntrdarray)) continue; // bug fix | |
294 | trdarray->AddPoint(ntrdarray++,&p); | |
295 | fCutpoi->Fill(15); | |
296 | fModpop->Fill(iModule,layer); | |
297 | ntrd++; | |
298 | } | |
299 | //printf("event %4d track %4d volid %4d layer %4d module %4d %10.3f %10.3f %10.3f\n", | |
300 | // (int)entry,itrack,volid,layer,iModule,p.GetX(),p.GetY(),atan2(p.GetY(),p.GetX())*180/3.1416); | |
301 | tmpArray->AddPoint(current,&p); | |
302 | current++; | |
303 | fCutpoi->Fill(16); | |
304 | } | |
305 | fCuttra->Fill(8); | |
306 | if (ntpc < 100) continue; | |
307 | fCuttra->Fill(9); | |
308 | if (ntrd < 4) continue; | |
309 | fCuttra->Fill(10); | |
310 | if (fArray) delete fArray; | |
311 | fArray = new AliTrackPointArray(current); | |
312 | for (Int_t ipoint=0; ipoint<current; ipoint++){ | |
313 | AliTrackPoint p; | |
314 | tmpArray->GetPoint(p, ipoint); | |
315 | fArray->AddPoint(ipoint,&p); | |
316 | } | |
317 | fNAcceptedTracks++; | |
318 | fTree->Fill(); | |
319 | } | |
320 | delete trdarray; | |
321 | fNevents++; | |
322 | return kTRUE; | |
323 | } | |
324 | //============================================================================= | |
325 | void PointMaker::SlaveTerminate() | |
326 | { | |
327 | // The SlaveTerminate() function is called after all entries or objects | |
328 | // have been processed. When running with PROOF SlaveTerminate() is called | |
329 | // on each slave server. | |
330 | ||
331 | // Add the histograms to the output on each slave server | |
332 | AliInfo(Form("\nNumber of tracks:\tprocessed\t%d\taccepted\t%d", fNtracks, fNAcceptedTracks)); | |
333 | } | |
334 | //============================================================================= | |
335 | void PointMaker::Terminate() | |
336 | { | |
337 | // The Terminate() function is the last function to be called during | |
338 | // a query. It always runs on the client, it can be used to present | |
339 | // the results graphically or save the results to file. | |
340 | ||
341 | // fTree = dynamic_cast<TTree*>(fOutput->FindObject("spTree")); | |
342 | ||
343 | // TFile* file = TFile::Open("AliTrackPoints.root", "RECREATE"); | |
344 | fFile->cd(); | |
345 | fTree->Write(); | |
346 | fCuttra->Write(); | |
347 | fCutpoi->Write(); | |
348 | fModpop->Write(); | |
349 | fModpop->Draw(); | |
350 | } | |
351 | //============================================================================= |