]>
Commit | Line | Data |
---|---|---|
f27a7e81 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-2009, 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 | ///////////////////////////////////////////////////////////// | |
17 | // | |
18 | // AliAnalysisTask to extract from ESD tracks the AliTrackPointArrays | |
19 | // with ITS points for selected tracks. This are the input data for alignment | |
20 | // | |
367c6d1f | 21 | // Author: A.Dainese, andrea.dainese@pd.infn.it |
f27a7e81 | 22 | ///////////////////////////////////////////////////////////// |
23 | ||
24 | #include <TTree.h> | |
eef875ec | 25 | #include <TFile.h> |
26 | #include <TSystem.h> | |
f27a7e81 | 27 | #include <TChain.h> |
28 | #include <TNtuple.h> | |
29 | #include <TList.h> | |
30 | #include <TH1F.h> | |
31 | #include <TH2F.h> | |
cf53e6db | 32 | #include <TMap.h> |
f27a7e81 | 33 | #include <TVector3.h> |
34 | #include <TGeoManager.h> | |
28647bf6 | 35 | #include <TRandom.h> |
f27a7e81 | 36 | |
37 | #include "AliLog.h" | |
eef875ec | 38 | #include "AliCDBEntry.h" |
f27a7e81 | 39 | #include "AliGeomManager.h" |
a043746c | 40 | #include "AliITSReconstructor.h" |
eef875ec | 41 | #include "AliITSAlignMille2Module.h" |
f27a7e81 | 42 | #include "AliITSgeomTGeo.h" |
43 | #include "AliTrackPointArray.h" | |
44 | #include "AliESDInputHandler.h" | |
45 | #include "AliESDVertex.h" | |
46 | #include "AliESDtrack.h" | |
47 | #include "AliESDEvent.h" | |
48 | #include "AliESDfriend.h" | |
b900a060 | 49 | #include "AliAnalysisTaskSE.h" |
f27a7e81 | 50 | #include "AliAnalysisManager.h" |
51 | #include "AliAlignmentDataFilterITS.h" | |
52 | ||
53 | ||
54 | ClassImp(AliAlignmentDataFilterITS) | |
55 | ||
56 | ||
b900a060 | 57 | //________________________________________________________________________ |
58 | AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(): | |
59 | AliAnalysisTaskSE(), | |
60 | fOnlySPDFO(kFALSE), | |
28647bf6 | 61 | fDownsamplelowpt(kFALSE), |
b900a060 | 62 | fGeometryFileName("geometry.root"), |
63 | fITSRecoParam(0), | |
64 | fESD(0), | |
65 | fListOfHistos(0), | |
66 | fspTree(0), | |
67 | fHistNevents(0), | |
68 | fHistNpoints(0), | |
69 | fHistPt(0), | |
70 | fHistLayer0(0), | |
71 | fHistLayer1(0), | |
72 | fHistLayer2(0), | |
73 | fHistLayer3(0), | |
74 | fHistLayer4(0), | |
75 | fHistLayer5(0), | |
76 | fntExtra(0), | |
77 | fntCosmicMatching(0) | |
78 | { | |
79 | // Constructor | |
80 | } | |
81 | ||
f27a7e81 | 82 | //________________________________________________________________________ |
83 | AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name): | |
b900a060 | 84 | AliAnalysisTaskSE(name), |
367c6d1f | 85 | fOnlySPDFO(kFALSE), |
28647bf6 | 86 | fDownsamplelowpt(kFALSE), |
367c6d1f | 87 | fGeometryFileName("geometry.root"), |
88 | fITSRecoParam(0), | |
f27a7e81 | 89 | fESD(0), |
f27a7e81 | 90 | fListOfHistos(0), |
91 | fspTree(0), | |
480207c6 | 92 | fHistNevents(0), |
f27a7e81 | 93 | fHistNpoints(0), |
94 | fHistPt(0), | |
95 | fHistLayer0(0), | |
96 | fHistLayer1(0), | |
97 | fHistLayer2(0), | |
98 | fHistLayer3(0), | |
99 | fHistLayer4(0), | |
100 | fHistLayer5(0), | |
101 | fntExtra(0), | |
102 | fntCosmicMatching(0) | |
103 | { | |
104 | // Constructor | |
105 | ||
b900a060 | 106 | // Define output slots here |
107 | ||
108 | // Output slot #1 writes into a TTree | |
109 | DefineOutput(1,TTree::Class()); //My private output | |
110 | // Output slot #2 writes into a TList | |
111 | DefineOutput(2,TList::Class()); //My private output | |
f27a7e81 | 112 | } |
113 | ||
114 | //________________________________________________________________________ | |
115 | AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS() | |
116 | { | |
117 | // Destructor | |
118 | if (fListOfHistos) { | |
119 | delete fListOfHistos; | |
120 | fListOfHistos = 0; | |
121 | } | |
122 | if (fspTree) { | |
123 | delete fspTree; | |
124 | fspTree = 0; | |
125 | } | |
480207c6 | 126 | if (fHistNevents) { |
127 | delete fHistNevents; | |
128 | fHistNevents = 0; | |
129 | } | |
a043746c | 130 | if (fHistNpoints) { |
131 | delete fHistNpoints; | |
132 | fHistNpoints = 0; | |
133 | } | |
f27a7e81 | 134 | if (fHistPt) { |
135 | delete fHistPt; | |
136 | fHistPt = 0; | |
137 | } | |
138 | if (fHistLayer0) { | |
139 | delete fHistLayer0; | |
140 | fHistLayer0 = 0; | |
141 | } | |
142 | if (fHistLayer1) { | |
143 | delete fHistLayer1; | |
144 | fHistLayer1 = 0; | |
145 | } | |
146 | if (fHistLayer2) { | |
147 | delete fHistLayer2; | |
148 | fHistLayer2 = 0; | |
149 | } | |
150 | if (fHistLayer3) { | |
151 | delete fHistLayer3; | |
152 | fHistLayer3 = 0; | |
153 | } | |
154 | if (fHistLayer4) { | |
155 | delete fHistLayer4; | |
156 | fHistLayer4 = 0; | |
157 | } | |
158 | if (fHistLayer5) { | |
159 | delete fHistLayer5; | |
160 | fHistLayer5 = 0; | |
161 | } | |
162 | if (fntExtra) { | |
163 | delete fntExtra; | |
164 | fntExtra = 0; | |
165 | } | |
166 | if (fntCosmicMatching) { | |
167 | delete fntCosmicMatching; | |
168 | fntCosmicMatching = 0; | |
169 | } | |
170 | } | |
a043746c | 171 | |
f27a7e81 | 172 | |
173 | //________________________________________________________________________ | |
b900a060 | 174 | void AliAlignmentDataFilterITS::UserCreateOutputObjects() |
f27a7e81 | 175 | { |
176 | // Create the output container | |
177 | // | |
c1605e49 | 178 | |
179 | // load the geometry | |
180 | if(!gGeoManager) { | |
181 | printf("AliAlignmentDataFilterITS::CreateOutputObjects(): loading geometry from %s\n",fGeometryFileName.Data()); | |
182 | AliGeomManager::LoadGeometry(fGeometryFileName.Data()); | |
183 | if(!gGeoManager) { | |
184 | printf("AliAlignmentDataFilterITS::CreateOutputObjects(): no geometry loaded \n"); | |
185 | return; | |
186 | } | |
187 | } | |
f27a7e81 | 188 | |
189 | // Several histograms are more conveniently managed in a TList | |
190 | fListOfHistos = new TList(); | |
191 | fListOfHistos->SetOwner(); | |
a043746c | 192 | |
480207c6 | 193 | fHistNevents = new TH1F("fHistNevents", "Number of processed events; N events; bin",5,-0.5,4.5); |
194 | fHistNevents->Sumw2(); | |
195 | fHistNevents->SetMinimum(0); | |
196 | fListOfHistos->Add(fHistNevents); | |
f27a7e81 | 197 | |
198 | fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5); | |
199 | fHistNpoints->Sumw2(); | |
200 | fHistNpoints->SetMinimum(0); | |
201 | fListOfHistos->Add(fHistNpoints); | |
202 | ||
203 | fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50); | |
204 | fHistPt->Sumw2(); | |
205 | fHistPt->SetMinimum(0); | |
206 | fListOfHistos->Add(fHistPt); | |
207 | ||
208 | ||
209 | Float_t zmax=14.; | |
210 | Int_t nbinsphi=20,nbinsz=4; | |
211 | fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); | |
212 | fListOfHistos->Add(fHistLayer0); | |
213 | zmax=14.; | |
214 | nbinsphi=40;nbinsz=4; | |
215 | fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); | |
216 | fListOfHistos->Add(fHistLayer1); | |
217 | zmax=22.; | |
218 | nbinsphi=14;nbinsz=6; | |
219 | fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); | |
220 | fListOfHistos->Add(fHistLayer2); | |
221 | zmax=29.5; | |
222 | nbinsphi=22;nbinsz=8; | |
223 | fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); | |
224 | fListOfHistos->Add(fHistLayer3); | |
225 | zmax=45.; | |
226 | nbinsphi=34;nbinsz=23; | |
227 | fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); | |
228 | fListOfHistos->Add(fHistLayer4); | |
229 | zmax=51.; | |
230 | nbinsphi=38;nbinsz=26; | |
231 | fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax); | |
232 | fListOfHistos->Add(fHistLayer5); | |
233 | ||
234 | ||
cdf8a7eb | 235 | fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu:pt"); |
f27a7e81 | 236 | fListOfHistos->Add(fntExtra); |
237 | ||
238 | fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu"); | |
239 | fListOfHistos->Add(fntCosmicMatching); | |
240 | ||
3acc122a | 241 | fspTree = new TTree("spTree","Tree with ITS track points"); |
eef875ec | 242 | AliTrackPointArray *array = 0; |
44f17b29 | 243 | AliESDVertex *vertex = 0; |
eef875ec | 244 | Float_t curv=0,curverr=0,runNumber=0; |
245 | TObjString *itsaligndata = 0; | |
246 | TObjString *itscalibrespsdd = 0; | |
f27a7e81 | 247 | fspTree->Branch("SP","AliTrackPointArray",&array); |
44f17b29 | 248 | fspTree->Branch("vertex","AliESDVertex",&vertex); |
f27a7e81 | 249 | fspTree->Branch("curv",&curv); |
250 | fspTree->Branch("curverr",&curverr); | |
ab6c74ff | 251 | fspTree->Branch("run",&runNumber); |
252 | fspTree->Branch("ITSAlignData",&itsaligndata); | |
3b9267a6 | 253 | fspTree->Branch("ITSCalibRespSDD",&itscalibrespsdd); |
f27a7e81 | 254 | |
255 | return; | |
256 | } | |
257 | ||
258 | //________________________________________________________________________ | |
b900a060 | 259 | void AliAlignmentDataFilterITS::UserExec(Option_t */*option*/) |
f27a7e81 | 260 | { |
261 | // Execute analysis for current event: | |
262 | // write ITS AliTrackPoints for selected tracks to fspTree | |
a043746c | 263 | |
c1605e49 | 264 | // check the geometry |
265 | if(!gGeoManager) { | |
266 | printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n"); | |
267 | return; | |
367c6d1f | 268 | } |
269 | ||
270 | // check if we have AliITSRecoParam | |
a043746c | 271 | if(!GetRecoParam()) { |
367c6d1f | 272 | if(!fITSRecoParam) { |
273 | printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n"); | |
274 | return; | |
275 | } | |
f27a7e81 | 276 | } |
277 | ||
b900a060 | 278 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); |
f27a7e81 | 279 | if(!fESD) { |
280 | printf("AliAlignmentDataFilterITS::Exec(): no ESD \n"); | |
281 | return; | |
282 | } | |
f27a7e81 | 283 | |
28647bf6 | 284 | //AliESDfriend *esdfriend = (AliESDfriend*)(fESD->FindListObject("AliESDfriend")); |
285 | ||
286 | //if(!esdfriend) printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n"); | |
287 | //fESD->SetESDfriend(esdfriend); | |
288 | ||
480207c6 | 289 | // Post the data for slot 0 |
290 | fHistNevents->Fill(0); | |
291 | ||
cf53e6db | 292 | |
293 | // write field value to spTree UserInfo | |
294 | if(!((fspTree->GetUserInfo())->FindObject("BzkGauss"))) { | |
295 | Double_t bz=fESD->GetMagneticField(); | |
296 | TString bzString; bzString+=bz; | |
297 | TObjString *bzObjString = new TObjString(bzString); | |
298 | TList *bzList = new TList(); | |
299 | bzList->SetOwner(1); | |
300 | bzList->SetName("BzkGauss"); | |
301 | bzList->Add(bzObjString); | |
302 | fspTree->GetUserInfo()->Add(bzList); | |
303 | } | |
304 | ||
b900a060 | 305 | // write OCDB info to spTree UserInfo |
306 | if(!((fspTree->GetUserInfo())->FindObject("cdbList"))) { | |
307 | TTree* tree = dynamic_cast<TTree*> (GetInputData(0)); | |
308 | if(!tree) { | |
309 | printf("ERROR: Could not read chain from input slot 0\n"); | |
310 | } else { | |
311 | // Get the OCDB path and the list of OCDB objects used for reco | |
312 | TMap *cdbMap = (TMap*)(tree->GetTree()->GetUserInfo())->FindObject("cdbMap"); | |
313 | TList *cdbList = (TList*)(tree->GetTree()->GetUserInfo())->FindObject("cdbList"); | |
314 | ||
315 | //cdbList->Print(); | |
316 | // write the list to the user info of the output tree | |
317 | if(!fspTree) { | |
318 | printf("ERROR: fspTree does not exist\n"); | |
319 | } else { | |
320 | TMap *cdbMapCopy = new TMap(cdbMap->GetEntries()); | |
321 | cdbMapCopy->SetOwner(1); | |
322 | cdbMapCopy->SetName("cdbMap"); | |
323 | TIter iter1(cdbMap->GetTable()); | |
324 | ||
325 | TPair* pair = 0; | |
326 | while((pair = dynamic_cast<TPair*> (iter1.Next()))){ | |
327 | TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key()); | |
328 | TObjString* valStr = dynamic_cast<TObjString*> (pair->Value()); | |
b6292968 | 329 | if(keyStr && valStr) cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName())); |
b900a060 | 330 | } |
331 | ||
332 | TList *cdbListCopy = new TList(); | |
333 | cdbListCopy->SetOwner(1); | |
334 | cdbListCopy->SetName("cdbList"); | |
335 | ||
336 | TIter iter2(cdbList); | |
337 | ||
338 | TObjString* cdbEntry=0; | |
339 | while((cdbEntry =(TObjString*)(iter2.Next()))) { | |
340 | cdbListCopy->Add(new TObjString(*cdbEntry)); | |
341 | } | |
342 | cdbListCopy->Print(); | |
343 | ||
344 | ||
345 | fspTree->GetUserInfo()->Add(cdbMapCopy); | |
346 | fspTree->GetUserInfo()->Add(cdbListCopy); | |
347 | } | |
348 | } | |
349 | } | |
350 | ||
351 | ||
cf53e6db | 352 | |
f27a7e81 | 353 | // Process event as Cosmic or Collision |
44f17b29 | 354 | if(fESD->GetEventSpecie()<=1) { |
355 | printf("AliAlignmentDataFilterITS::Exec(): event specie not set !\n"); | |
356 | if(GetRecoParam()->GetAlignFilterCosmics()) { | |
357 | FilterCosmic(fESD); | |
358 | } else { | |
359 | FilterCollision(fESD); | |
360 | } | |
361 | } else if(fESD->GetEventSpecie()==8) { | |
f27a7e81 | 362 | FilterCosmic(fESD); |
363 | } else { | |
364 | FilterCollision(fESD); | |
365 | } | |
366 | ||
b900a060 | 367 | PostData(2,fListOfHistos); |
480207c6 | 368 | |
f27a7e81 | 369 | return; |
370 | } | |
371 | ||
372 | //________________________________________________________________________ | |
373 | void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd) | |
374 | { | |
375 | // Extract ITS AliTrackPoints for Cosmics (check angular matching | |
376 | // of top and bottom track, merge the two tracks, if requested) | |
377 | // | |
378 | ||
379 | // Set branch addresses for space points tree | |
380 | AliTrackPointArray *arrayForTree=0; | |
44f17b29 | 381 | AliESDVertex *vertexForTree=0; |
ab6c74ff | 382 | Float_t curv,curverr,runNumber; |
383 | TObjString *itsaligndata=0; | |
3b9267a6 | 384 | TObjString *itscalibrespsdd = 0; |
f27a7e81 | 385 | fspTree->SetBranchAddress("SP",&arrayForTree); |
44f17b29 | 386 | fspTree->SetBranchAddress("vertex",&vertexForTree); |
f27a7e81 | 387 | fspTree->SetBranchAddress("curv",&curv); |
388 | fspTree->SetBranchAddress("curverr",&curverr); | |
ab6c74ff | 389 | fspTree->SetBranchAddress("run",&runNumber); |
390 | fspTree->SetBranchAddress("ITSAlignData",&itsaligndata); | |
3b9267a6 | 391 | fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd); |
ab6c74ff | 392 | |
393 | ||
394 | runNumber = (Float_t)esd->GetRunNumber(); | |
b900a060 | 395 | Int_t uid=10000+esd->GetEventNumberInFile(); |
ab6c74ff | 396 | |
397 | TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0)); | |
b6292968 | 398 | if(!esdTree) return; |
ab6c74ff | 399 | // Get the list of OCDB objects used for reco |
400 | TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList"); | |
401 | TIter iter2(cdbList); | |
402 | TObjString* cdbEntry=0; | |
403 | TString cdbEntryString; | |
404 | while((cdbEntry =(TObjString*)(iter2.Next()))) { | |
ca6660ea | 405 | cdbEntryString = cdbEntry->GetString(); |
406 | if(cdbEntryString.Contains("ITS/Align/Data")) { | |
407 | itsaligndata = new TObjString(*cdbEntry); | |
408 | itsaligndata->SetString(itsaligndata->GetString()); | |
409 | } | |
410 | if(cdbEntryString.Contains("ITS/Calib/RespSDD")) { | |
411 | itscalibrespsdd = new TObjString(*cdbEntry); | |
412 | itscalibrespsdd->SetString(itscalibrespsdd->GetString()); | |
413 | } | |
ab6c74ff | 414 | } |
415 | ||
f27a7e81 | 416 | |
b900a060 | 417 | TString triggeredClass = esd->GetFiredTriggerClasses(); |
367c6d1f | 418 | if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return; |
419 | ||
f27a7e81 | 420 | |
421 | Int_t ntracks = esd->GetNumberOfTracks(); | |
422 | if(ntracks<2) return; | |
423 | ||
424 | if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return; | |
425 | ||
426 | Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos); | |
427 | ||
428 | Int_t *goodtracksArray = new Int_t[ntracks]; | |
429 | Float_t *phiArray = new Float_t[ntracks]; | |
430 | Float_t *thetaArray = new Float_t[ntracks]; | |
431 | Int_t *nclsArray = new Int_t[ntracks]; | |
432 | Int_t ngt=0; | |
433 | Int_t itrack=0; | |
434 | for (itrack=0; itrack < ntracks; itrack++) { | |
435 | AliESDtrack *track = esd->GetTrack(itrack); | |
436 | if (!track) continue; | |
437 | ||
438 | ||
367c6d1f | 439 | if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue; |
f27a7e81 | 440 | |
a043746c | 441 | if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue; |
442 | if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue; | |
f27a7e81 | 443 | |
444 | Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp()); | |
445 | Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl()); | |
446 | ||
367c6d1f | 447 | if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() || |
448 | track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue; | |
f27a7e81 | 449 | |
450 | goodtracksArray[ngt] = itrack; | |
451 | phiArray[ngt] = phi; | |
452 | thetaArray[ngt] = theta; | |
453 | nclsArray[ngt] = track->GetNcls(0); | |
454 | ngt++; | |
455 | } | |
456 | ||
457 | if(ngt<2) { | |
458 | delete [] goodtracksArray; goodtracksArray=0; | |
459 | delete [] phiArray; phiArray=0; | |
460 | delete [] thetaArray; thetaArray=0; | |
461 | delete [] nclsArray; nclsArray=0; | |
462 | return; | |
463 | } | |
464 | ||
465 | // check matching of the two tracks from the muon | |
466 | Float_t min = 10000000.; | |
467 | Int_t maxCls = 0; | |
468 | Int_t good1 = -1, good2 = -1; | |
469 | for(Int_t itr1=0; itr1<ngt-1; itr1++) { | |
470 | for(Int_t itr2=itr1+1; itr2<ngt; itr2++) { | |
471 | Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]); | |
367c6d1f | 472 | if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue; |
f27a7e81 | 473 | Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi()); |
367c6d1f | 474 | if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue; |
f27a7e81 | 475 | if(nclsArray[itr1]+nclsArray[itr2] > maxCls) { |
476 | maxCls = nclsArray[itr1]+nclsArray[itr2]; | |
477 | min = deltatheta+deltaphi; | |
478 | good1 = goodtracksArray[itr1]; | |
479 | good2 = goodtracksArray[itr2]; | |
480 | } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) { | |
481 | if(deltatheta+deltaphi < min) { | |
482 | min = deltatheta+deltaphi; | |
483 | good1 = goodtracksArray[itr1]; | |
484 | good2 = goodtracksArray[itr2]; | |
485 | } | |
486 | } | |
487 | } | |
488 | } | |
489 | ||
490 | delete [] goodtracksArray; goodtracksArray=0; | |
491 | delete [] phiArray; phiArray=0; | |
492 | delete [] thetaArray; thetaArray=0; | |
493 | delete [] nclsArray; nclsArray=0; | |
494 | ||
495 | if(good1<0) return; | |
496 | AliDebug(2,"ok track matching"); | |
497 | ||
498 | // track1 will be the inward track (top) | |
499 | // track2 the outward (bottom) | |
500 | AliESDtrack *track1=0; | |
501 | AliESDtrack *track2=0; | |
502 | AliESDtrack *track = esd->GetTrack(good1); | |
503 | if(track->Py()>0) { | |
504 | track1 = esd->GetTrack(good1); | |
505 | track2 = esd->GetTrack(good2); | |
506 | } else { | |
507 | track1 = esd->GetTrack(good2); | |
508 | track2 = esd->GetTrack(good1); | |
509 | } | |
510 | ||
511 | AliTrackPoint point; | |
512 | const AliTrackPointArray *array=0; | |
513 | Int_t ipt,volId,modId,layerId,lay,lad,det; | |
514 | Int_t jpt=0; | |
515 | Bool_t layerOK[6][2]; | |
516 | Int_t nclsTrk[2]={0,0}; | |
517 | ||
518 | for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE; | |
519 | ||
520 | for(itrack=0; itrack<2; itrack++) { | |
521 | if(itrack==0) { | |
522 | track = track1; | |
523 | } else { | |
524 | track = track2; | |
525 | } | |
526 | array = track->GetTrackPointArray(); | |
527 | if(!array) { | |
528 | AliWarning("No tracks points avaialble"); | |
529 | continue; | |
530 | } | |
531 | for(ipt=0; ipt<array->GetNPoints(); ipt++) { | |
532 | array->GetPoint(point,ipt); | |
533 | volId = point.GetVolumeID(); | |
9f2b55aa | 534 | if(volId<=0) continue; |
f27a7e81 | 535 | layerId = AliGeomManager::VolUIDToLayer(volId,modId); |
9f2b55aa | 536 | AliDebug(2,Form("%d %d %d %f\n",ipt,layerId-1,volId,TMath::Sqrt(point.GetX()*point.GetX()+point.GetY()*point.GetY()))); |
f27a7e81 | 537 | if(point.IsExtra() && |
a043746c | 538 | (GetRecoParam()->GetAlignFilterSkipExtra())) continue; |
9f2b55aa | 539 | if(layerId<1 || layerId>6) continue; |
367c6d1f | 540 | if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue; |
f27a7e81 | 541 | // check minAngleWrtITSModulePlanes |
542 | Double_t p[3]; track->GetDirection(p); | |
543 | TVector3 pvec(p[0],p[1],p[2]); | |
544 | Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot); | |
545 | TVector3 normvec(rot[1],rot[4],rot[7]); | |
546 | Double_t angle = pvec.Angle(normvec); | |
547 | if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle; | |
548 | angle = 0.5*TMath::Pi()-angle; | |
367c6d1f | 549 | if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue; |
f27a7e81 | 550 | layerOK[layerId-1][itrack]=kTRUE; |
551 | jpt++; | |
552 | nclsTrk[itrack]++; | |
553 | } | |
554 | } | |
555 | AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1])); | |
556 | ||
557 | // read ITS cluster maps | |
558 | Int_t map1[6],map2[6]; | |
559 | for(Int_t ilay=0;ilay<6;ilay++) { | |
560 | map1[ilay]=0; map2[ilay]=0; | |
561 | if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1; | |
562 | if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1; | |
563 | } | |
564 | AliDebug(2,Form("ITS map 1: %d %d %d %d %d %d pt %f\n",map1[0],map1[1],map1[2],map1[3],map1[4],map1[5],track1->Pt())); | |
565 | AliDebug(2,Form("ITS map 2: %d %d %d %d %d %d pt %f\n",map2[0],map2[1],map2[2],map2[3],map2[4],map2[5],track2->Pt())); | |
566 | Int_t idx1[12],idx2[12]; | |
567 | track1->GetITSclusters(idx1); | |
568 | track2->GetITSclusters(idx2); | |
569 | AliDebug(2,Form("cls idx 1 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx1[0],idx1[1],idx1[2],idx1[3],idx1[4],idx1[5],idx1[6],idx1[7],idx1[8],idx1[9],idx1[10],idx1[11])); | |
570 | AliDebug(2,Form("cls idx 2 %d %d %d %d %d %d %d %d %d %d %d %d\n",idx2[0],idx2[1],idx2[2],idx2[3],idx2[4],idx2[5],idx2[6],idx2[7],idx2[8],idx2[9],idx2[10],idx2[11])); | |
571 | ||
572 | ||
367c6d1f | 573 | if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return; |
f27a7e81 | 574 | AliDebug(2,Form(" Total points %d, accepted\n",jpt)); |
575 | fHistNpoints->Fill(jpt); | |
576 | fHistPt->Fill(0.5*(track1->Pt()+track2->Pt())); | |
577 | ||
578 | Float_t d0z0mu[2]; | |
579 | track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu); | |
580 | //printf("d0mu %f z0mu %f\n",d0z0mu[0],d0z0mu[1]); | |
581 | ||
44f17b29 | 582 | vertexForTree = new AliESDVertex(*(esd->GetPrimaryVertexSPD())); |
b6292968 | 583 | vertexForTree->SetID(0); |
44f17b29 | 584 | |
f27a7e81 | 585 | Float_t dzOverlap[2]; |
586 | Float_t curvArray[2],curverrArray[2]; | |
587 | Double_t globExtra[3],locExtra[3]; | |
cdf8a7eb | 588 | if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) { |
f27a7e81 | 589 | arrayForTree = new AliTrackPointArray(jpt); |
cdf8a7eb | 590 | arrayForTree->SetUniqueID(uid); |
591 | } | |
f27a7e81 | 592 | jpt=0; |
593 | for(itrack=0; itrack<2; itrack++) { | |
594 | if(itrack==0) { | |
595 | track = track1; | |
596 | } else { | |
597 | track = track2; | |
598 | } | |
599 | curvArray[itrack] = track->GetC(esd->GetMagneticField()); | |
600 | curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt(); | |
601 | ||
a043746c | 602 | if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) { |
f27a7e81 | 603 | jpt=0; |
604 | arrayForTree = new AliTrackPointArray(nclsTrk[itrack]); | |
cdf8a7eb | 605 | arrayForTree->SetUniqueID(uid); |
f27a7e81 | 606 | } |
607 | array = track->GetTrackPointArray(); | |
608 | for(ipt=0; ipt<array->GetNPoints(); ipt++) { | |
609 | array->GetPoint(point,ipt); | |
610 | volId = point.GetVolumeID(); | |
9f2b55aa | 611 | if(volId<=0) continue; |
f27a7e81 | 612 | layerId = AliGeomManager::VolUIDToLayer(volId,modId); |
9f2b55aa | 613 | if(layerId<1 || layerId>6 || !layerOK[layerId-1][itrack]) continue; |
f27a7e81 | 614 | arrayForTree->AddPoint(jpt,&point); |
615 | jpt++; | |
616 | switch(layerId) { | |
617 | case 1: | |
618 | fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
619 | break; | |
620 | case 2: | |
621 | fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
622 | break; | |
623 | case 3: | |
624 | fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
625 | break; | |
626 | case 4: | |
627 | fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
628 | break; | |
629 | case 5: | |
630 | fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
631 | break; | |
632 | case 6: | |
633 | fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
634 | break; | |
635 | } | |
b900a060 | 636 | // Post the data for slot 2 |
637 | if(jpt==1) PostData(2,fListOfHistos); // only if this is the first points | |
f27a7e81 | 638 | if(!point.IsExtra() || |
a043746c | 639 | !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue; |
f27a7e81 | 640 | nclsTrk[itrack]--; |
641 | for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll); | |
642 | AliITSgeomTGeo::GetModuleId(modId,lay,lad,det); | |
643 | globExtra[0]=point.GetX(); | |
644 | globExtra[1]=point.GetY(); | |
645 | globExtra[2]=point.GetZ(); | |
646 | AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra); | |
647 | //printf("%d %d %d %d %d %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]); | |
648 | track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap); | |
649 | AliTrackPoint pointT; | |
650 | Float_t radius,radiusT,phiv,phivT,thetav,thetavT; | |
651 | for(Int_t lll=0;lll<ipt;lll++) { | |
652 | array->GetPoint(pointT,lll); | |
9f2b55aa | 653 | if(pointT.GetVolumeID()<=0) continue; |
f27a7e81 | 654 | Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId); |
655 | if(layerIdT!=layerId) continue; | |
656 | radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1])); | |
657 | radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1])); | |
658 | phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]); | |
659 | phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]); | |
660 | if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue; | |
661 | thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2])); | |
662 | thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2])); | |
663 | dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius)); | |
664 | if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue; | |
665 | dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT))); | |
cdf8a7eb | 666 | fntExtra->Fill((Float_t)nclsTrk[itrack],(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0mu[0],d0z0mu[1],track->Pt()); |
f27a7e81 | 667 | } |
668 | } | |
669 | ||
a043746c | 670 | if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) { |
f27a7e81 | 671 | curv = curvArray[itrack]; |
672 | curverr = curverrArray[itrack]; | |
673 | fspTree->Fill(); | |
674 | } | |
675 | } | |
676 | ||
367c6d1f | 677 | if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) { |
3e33351b | 678 | curv = 0.5*(curvArray[0]-curvArray[1]); // the "-" is because the two tracks have opposite curvature! |
f27a7e81 | 679 | curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]); |
680 | fspTree->Fill(); | |
681 | } | |
b900a060 | 682 | PostData(1,fspTree); |
f27a7e81 | 683 | |
a043746c | 684 | if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return; |
f27a7e81 | 685 | // fill ntuple with track-to-track matching |
686 | Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty; | |
687 | Float_t d0[2],z0[2]; | |
688 | Float_t sigmad0[2],sigmaz0[2]; | |
689 | phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp()); | |
690 | thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl()); | |
691 | phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp()); | |
692 | thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl()); | |
693 | rotymu = TMath::ATan2(track1->Px(),track1->Pz()); | |
694 | rotyout = TMath::ATan2(track2->Px(),track2->Pz()); | |
695 | ||
696 | dphi = phimu - (phiout+TMath::Pi()); | |
697 | dtheta = thetamu - (TMath::Pi()-thetaout); | |
698 | if(rotymu>0) { | |
699 | droty = rotymu - (rotyout+TMath::Pi()); | |
700 | } else { | |
701 | droty = rotymu - (rotyout-TMath::Pi()); | |
702 | } | |
703 | ||
704 | Double_t alpha = TMath::ATan2(track1->Py(),track1->Px()); | |
705 | ||
706 | track1->Propagate(alpha,0.,esd->GetMagneticField()); | |
707 | track2->Propagate(alpha,0.,esd->GetMagneticField()); | |
708 | d0[0] = track1->GetY(); | |
709 | z0[0] = track1->GetZ(); | |
710 | d0[1] = track2->GetY(); | |
711 | z0[1] = track2->GetZ(); | |
712 | Float_t dxy = -(d0[0]-d0[1]); | |
713 | Float_t dz = z0[0]-z0[1]; | |
714 | sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2()); | |
715 | sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2()); | |
716 | sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2()); | |
717 | sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2()); | |
718 | /* | |
719 | Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3]; | |
720 | track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0); | |
721 | track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1); | |
722 | track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0); | |
723 | track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1); | |
724 | Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]); | |
725 | Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]); | |
726 | Float_t dxaty0 = x1aty0-x2aty0; | |
727 | */ | |
728 | fntCosmicMatching->Fill((Float_t)nclsTrk[0],(Float_t)nclsTrk[1],track1->Pt(),track2->Pt(),sigmad0[0],sigmad0[1],sigmaz0[0],sigmaz0[1],dxy,dz,phimu,thetamu,TMath::Abs(d0z0mu[0]),d0z0mu[1]); | |
729 | ||
730 | return; | |
731 | } | |
732 | ||
733 | //________________________________________________________________________ | |
734 | void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd) | |
735 | { | |
736 | // Extract ITS AliTrackPoints for Cosmics (check angular matching | |
737 | // of top and bottom track, merge the two tracks, if requested) | |
738 | // | |
739 | ||
740 | // Set branch addresses for space points tree | |
741 | AliTrackPointArray *arrayForTree=0; | |
44f17b29 | 742 | AliESDVertex *vertexForTree=0; |
ab6c74ff | 743 | Float_t curv,curverr,runNumber; |
744 | TObjString *itsaligndata=0; | |
3b9267a6 | 745 | TObjString *itscalibrespsdd = 0; |
f27a7e81 | 746 | fspTree->SetBranchAddress("SP",&arrayForTree); |
44f17b29 | 747 | fspTree->SetBranchAddress("vertex",&vertexForTree); |
f27a7e81 | 748 | fspTree->SetBranchAddress("curv",&curv); |
749 | fspTree->SetBranchAddress("curverr",&curverr); | |
ab6c74ff | 750 | fspTree->SetBranchAddress("run",&runNumber); |
751 | fspTree->SetBranchAddress("ITSAlignData",&itsaligndata); | |
3b9267a6 | 752 | fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd); |
ab6c74ff | 753 | |
754 | ||
755 | runNumber = (Float_t)esd->GetRunNumber(); | |
b900a060 | 756 | Int_t uid=20000+esd->GetEventNumberInFile(); |
ab6c74ff | 757 | |
758 | TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0)); | |
b6292968 | 759 | if(!esdTree) return; |
ab6c74ff | 760 | // Get the list of OCDB objects used for reco |
761 | TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList"); | |
762 | TIter iter2(cdbList); | |
763 | TObjString* cdbEntry=0; | |
764 | TString cdbEntryString; | |
765 | while((cdbEntry =(TObjString*)(iter2.Next()))) { | |
ca6660ea | 766 | cdbEntryString = cdbEntry->GetString(); |
767 | if(cdbEntryString.Contains("ITS/Align/Data")) { | |
768 | itsaligndata = new TObjString(*cdbEntry); | |
769 | itsaligndata->SetString(itsaligndata->GetString()); | |
770 | } | |
771 | if(cdbEntryString.Contains("ITS/Calib/RespSDD")) { | |
772 | itscalibrespsdd = new TObjString(*cdbEntry); | |
773 | itscalibrespsdd->SetString(itscalibrespsdd->GetString()); | |
774 | } | |
ab6c74ff | 775 | } |
f27a7e81 | 776 | |
777 | Int_t ntracks = esd->GetNumberOfTracks(); | |
778 | ||
779 | if(ntracks==0) return; | |
780 | ||
44f17b29 | 781 | const AliESDVertex *vertexTracks = esd->GetPrimaryVertexTracks(); |
9a4d2d47 | 782 | if(!vertexTracks) return; |
44f17b29 | 783 | if(vertexTracks->GetNContributors()<=0) return; |
f27a7e81 | 784 | |
44f17b29 | 785 | Double_t vtxpos[3]; vertexTracks->GetXYZ(vtxpos); |
f27a7e81 | 786 | |
787 | Int_t ncls=0; | |
788 | Double_t pt=-10000.; | |
789 | Double_t d0z0[2],covd0z0[3]; | |
790 | const AliTrackPointArray *array = 0; | |
791 | ||
792 | for (Int_t itrack=0; itrack < ntracks; itrack++) { | |
793 | AliESDtrack * track = esd->GetTrack(itrack); | |
794 | if (!track) continue; | |
795 | ||
28647bf6 | 796 | if(fDownsamplelowpt && TMath::Abs(esd->GetMagneticField())>0.01 && |
797 | track->Pt()<gRandom->Rndm()) continue; | |
798 | ||
367c6d1f | 799 | if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue; |
f27a7e81 | 800 | |
a043746c | 801 | if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue; |
802 | if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue; | |
f27a7e81 | 803 | |
367c6d1f | 804 | if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() || |
805 | track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue; | |
f27a7e81 | 806 | |
807 | pt = track->Pt(); | |
808 | ncls = track->GetNcls(0); | |
809 | Double_t maxd=10000.; | |
44f17b29 | 810 | track->PropagateToDCA(vertexTracks,esd->GetMagneticField(),maxd,d0z0,covd0z0); |
f27a7e81 | 811 | |
812 | // read ITS cluster map | |
813 | Int_t map[6]; | |
814 | for(Int_t ilay=0;ilay<6;ilay++) { | |
815 | map[ilay]=0; | |
816 | if(track->HasPointOnITSLayer(ilay)) map[ilay]=1; | |
817 | } | |
818 | AliDebug(2,Form("ITS map : %d %d %d %d %d %d pt %f\n",map[0],map[1],map[2],map[3],map[4],map[5],track->Pt())); | |
819 | Int_t idx[12]; | |
820 | track->GetITSclusters(idx); | |
821 | AliDebug(2,Form("cls idx %d %d %d %d %d %d %d %d %d %d %d %d\n",idx[0],idx[1],idx[2],idx[3],idx[4],idx[5],idx[6],idx[7],idx[8],idx[9],idx[10],idx[11])); | |
822 | ||
823 | ||
824 | AliTrackPoint point; | |
825 | Int_t ipt,volId,modId,layerId,lay,lad,det; | |
826 | Int_t jpt=0; | |
827 | Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE; | |
828 | ||
829 | array = track->GetTrackPointArray(); | |
28647bf6 | 830 | if(!array) {printf("no track points\n"); continue;} |
f27a7e81 | 831 | for(ipt=0; ipt<array->GetNPoints(); ipt++) { |
832 | array->GetPoint(point,ipt); | |
833 | volId = point.GetVolumeID(); | |
9f2b55aa | 834 | if(volId<=0) continue; |
f27a7e81 | 835 | layerId = AliGeomManager::VolUIDToLayer(volId,modId); |
836 | if(layerId<1 || layerId>6) continue; | |
837 | if(point.IsExtra() && | |
a043746c | 838 | (GetRecoParam()->GetAlignFilterSkipExtra())) continue; |
f27a7e81 | 839 | layerOK[layerId-1]=kTRUE; |
840 | jpt++; | |
841 | } | |
842 | ||
367c6d1f | 843 | if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue; |
f27a7e81 | 844 | |
845 | fHistNpoints->Fill(jpt); | |
846 | fHistPt->Fill(pt); | |
b900a060 | 847 | PostData(2,fListOfHistos); |
f27a7e81 | 848 | |
849 | Float_t dzOverlap[2]; | |
850 | Double_t globExtra[3],locExtra[3]; | |
851 | arrayForTree = new AliTrackPointArray(jpt); | |
cdf8a7eb | 852 | arrayForTree->SetUniqueID(uid); |
f27a7e81 | 853 | jpt=0; |
854 | array = track->GetTrackPointArray(); | |
855 | if(!array) continue; | |
856 | for(ipt=0; ipt<array->GetNPoints(); ipt++) { | |
857 | array->GetPoint(point,ipt); | |
858 | volId = point.GetVolumeID(); | |
859 | layerId = AliGeomManager::VolUIDToLayer(volId,modId); | |
860 | if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue; | |
44c3f5f6 | 861 | arrayForTree->AddPoint(jpt,&point); |
862 | switch(layerId) { | |
863 | case 1: | |
864 | fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
865 | break; | |
866 | case 2: | |
867 | fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
868 | break; | |
869 | case 3: | |
870 | fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
871 | break; | |
872 | case 4: | |
873 | fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
874 | break; | |
875 | case 5: | |
876 | fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
877 | break; | |
878 | case 6: | |
879 | fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ()); | |
880 | break; | |
881 | } | |
882 | jpt++; | |
f27a7e81 | 883 | if(!point.IsExtra() || |
a043746c | 884 | !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue; |
f27a7e81 | 885 | ncls--; |
886 | for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll); | |
887 | AliITSgeomTGeo::GetModuleId(modId,lay,lad,det); | |
888 | globExtra[0]=point.GetX(); | |
889 | globExtra[1]=point.GetY(); | |
890 | globExtra[2]=point.GetZ(); | |
891 | AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra); | |
892 | track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap); | |
893 | AliTrackPoint pointT; | |
894 | Float_t radius,radiusT,phiv,phivT,thetav,thetavT; | |
895 | for(Int_t lll=0;lll<ipt;lll++) { | |
896 | array->GetPoint(pointT,lll); | |
897 | Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId); | |
898 | if(layerIdT!=layerId) continue; | |
899 | radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1])); | |
900 | radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1])); | |
901 | phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]); | |
902 | phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]); | |
903 | if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue; | |
904 | thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2])); | |
905 | thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2])); | |
906 | dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius)); | |
907 | if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue; | |
908 | dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT))); | |
cdf8a7eb | 909 | fntExtra->Fill((Float_t)ncls,(Float_t)(layerId-1),lad,volId,TMath::ATan2(point.GetY(),point.GetX()),point.GetX(),point.GetY(),point.GetZ(),locExtra[0],locExtra[2],dzOverlap[0],dzOverlap[1],d0z0[0],d0z0[1],track->Pt()); |
f27a7e81 | 910 | } |
f27a7e81 | 911 | } |
912 | ||
913 | curv = track->GetC(esd->GetMagneticField()); | |
914 | curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt(); | |
915 | ||
44f17b29 | 916 | vertexForTree = new AliESDVertex(*vertexTracks); |
917 | if(vertexTracks->UsesTrack(track->GetID())) { | |
918 | vertexForTree->SetID(1); | |
919 | } else { | |
920 | vertexForTree->SetID(0); | |
921 | } | |
922 | ||
f27a7e81 | 923 | fspTree->Fill(); |
924 | ||
925 | } // end of tracks loop | |
926 | ||
b900a060 | 927 | PostData(1,fspTree); |
f27a7e81 | 928 | |
929 | return; | |
930 | } | |
931 | ||
932 | //________________________________________________________________________ | |
933 | void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/) | |
934 | { | |
935 | // Terminate analysis | |
936 | // | |
937 | AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n"); | |
938 | ||
81e84469 | 939 | fspTree = dynamic_cast<TTree*> (GetOutputData(1)); |
f27a7e81 | 940 | if (!fspTree) { |
941 | printf("ERROR: fspTree not available\n"); | |
942 | return; | |
943 | } | |
944 | ||
81e84469 | 945 | fListOfHistos = dynamic_cast<TList*> (GetOutputData(2)); |
f27a7e81 | 946 | if (!fListOfHistos) { |
947 | printf("ERROR: fListOfHistos not available\n"); | |
948 | return; | |
949 | } | |
950 | ||
a043746c | 951 | fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents")); |
f27a7e81 | 952 | fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints")); |
953 | fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt")); | |
954 | fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0")); | |
955 | fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1")); | |
956 | fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2")); | |
957 | fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3")); | |
958 | fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4")); | |
959 | fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5")); | |
960 | fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra")); | |
961 | fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching")); | |
962 | ||
963 | ||
964 | ||
965 | return; | |
966 | } | |
eef875ec | 967 | //------------------------------------------------------------------------------- |
968 | const AliITSRecoParam *AliAlignmentDataFilterITS::GetRecoParam() const | |
969 | { | |
970 | // | |
971 | // Return the ITSRecoParam object | |
972 | // | |
973 | if(AliITSReconstructor::GetRecoParam()) { | |
974 | return AliITSReconstructor::GetRecoParam(); | |
975 | } else if(fITSRecoParam) { | |
976 | return fITSRecoParam; | |
977 | } else return NULL; | |
978 | } | |
979 | //-------------------------------------------------------------------------------- | |
980 | Int_t AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom(Char_t *fin, | |
981 | Char_t *fout, | |
982 | Char_t *fmis, | |
983 | Char_t *fgeo, | |
984 | Bool_t prn) | |
985 | { | |
986 | // | |
987 | // Convert AliTrackPoints in fin, reconstructed with fmis, back | |
988 | // to ideal geometry | |
989 | // | |
990 | // M. Lunardon | |
991 | // | |
992 | ||
993 | ||
994 | TGeoHMatrix deltahm; | |
eef875ec | 995 | |
996 | // Load geometry | |
997 | if (gSystem->AccessPathName(fgeo)) { | |
998 | printf("couldn't find geometry file %s - skipping...\n",fmis); | |
999 | return -1; | |
1000 | } | |
1001 | ||
1002 | TFile *geofile=TFile::Open(fgeo); | |
1003 | TGeoManager *fgGeometry=NULL; | |
1004 | ||
1005 | fgGeometry=(TGeoManager*)geofile->Get("ALICE"); | |
1006 | ||
1007 | if (!fgGeometry) | |
1008 | fgGeometry=(TGeoManager*)geofile->Get("Geometry"); | |
1009 | ||
1010 | if (!fgGeometry) { | |
1011 | AliCDBEntry *entry = (AliCDBEntry*)geofile->Get("AliCDBEntry"); | |
1012 | if (entry) | |
1013 | fgGeometry = (TGeoManager*)entry->GetObject(); | |
1014 | } | |
1015 | ||
1016 | if (!fgGeometry) return -1; | |
1017 | AliGeomManager::SetGeometry(fgGeometry); | |
1018 | if(!AliGeomManager::GetGeometry()) return -1; | |
1019 | ||
1020 | ||
1021 | // open alignment file | |
1022 | if (gSystem->AccessPathName(fmis)) { | |
1023 | printf("couldn't open alignment file %s - skipping...\n",fmis); | |
1024 | return -2; | |
1025 | } | |
1026 | TFile *pref = TFile::Open(fmis); | |
1027 | if (!pref->IsOpen()) return -2; | |
1028 | ||
1029 | ||
1030 | /// apply alignment to ideal geometry | |
1031 | TClonesArray *prea=(TClonesArray*)pref->Get("ITSAlignObjs"); | |
1032 | if (!prea) { | |
1033 | if (pref->Get("AliCDBEntry")) | |
1034 | prea = (TClonesArray*) ((AliCDBEntry*)pref->Get("AliCDBEntry"))->GetObject(); | |
1035 | } | |
1036 | if (!prea) return -3; | |
1037 | Int_t nprea=prea->GetEntriesFast(); | |
1038 | printf("Array of input misalignments with %d entries\n",nprea); | |
1039 | AliGeomManager::ApplyAlignObjsToGeom(*prea); // apply all levels of objs | |
1040 | ||
1041 | AliTrackPointArray *tpain=NULL; | |
1042 | TFile *tpainfile=NULL; | |
1043 | TTree *treein=NULL; | |
1044 | AliTrackPoint point; | |
1045 | AliITSAlignMille2Module *m2[2200]; | |
1046 | for (Int_t i=0; i<2198; i++) | |
1047 | m2[i]=new AliITSAlignMille2Module(AliITSAlignMille2Module::GetVolumeIDFromIndex(i)); | |
1048 | ||
1049 | // open input file | |
1050 | if (gSystem->AccessPathName(fin)) { | |
1051 | printf("couldn't open file %s - skipping...\n",fin); | |
1052 | return -4; | |
1053 | } | |
1054 | tpainfile = TFile::Open(fin); | |
1055 | if (!tpainfile->IsOpen()) return -4; | |
1056 | ||
1057 | treein=(TTree*)tpainfile->Get("spTree"); | |
1058 | if (!treein) return -5; | |
1059 | Float_t curv,curverr,runNumber; | |
1060 | TObjString *itsaligndata=0; | |
1061 | TObjString *itscalibrespsdd = 0; | |
1062 | treein->SetBranchAddress("SP", &tpain); | |
1063 | treein->SetBranchAddress("curv", &curv); | |
1064 | treein->SetBranchAddress("curverr", &curverr); | |
1065 | treein->SetBranchAddress("run",&runNumber); | |
1066 | treein->SetBranchAddress("ITSAlignData",&itsaligndata); | |
1067 | treein->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd); | |
1068 | ||
1069 | int ntrks=treein->GetEntries(); | |
1070 | printf("Reading %d tracks from %s\n",ntrks,fin); | |
1071 | ||
1072 | ||
1073 | // open output file | |
1074 | TFile *pointsFile = TFile::Open(fout,"RECREATE"); | |
1075 | if (!pointsFile || !pointsFile->IsOpen()) { | |
1076 | printf("Can't open output file %s !",fout); | |
1077 | return -6; | |
1078 | } | |
1079 | AliTrackPointArray *array = new AliTrackPointArray(); | |
1080 | ||
1081 | // new! | |
1082 | TTree *treeout=(TTree*)treein->Clone("spTree"); | |
1083 | treeout->Reset(); | |
1084 | treeout->SetBranchAddress("SP", &array); | |
1085 | treeout->SetBranchAddress("curv", &curv); | |
1086 | treeout->SetBranchAddress("curverr", &curverr); | |
1087 | treeout->SetBranchAddress("run",&runNumber); | |
1088 | treeout->SetBranchAddress("ITSAlignData",&itsaligndata); | |
1089 | treeout->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd); | |
1090 | ||
1091 | // tracks main loop | |
1092 | for (Int_t it=0; it<ntrks; it++) { | |
1093 | if (!(it%5000) ) printf("...processing track n. %d\n",it); | |
1094 | ||
1095 | treein->GetEvent(it); | |
1096 | ||
1097 | ////////////////////////////// | |
1098 | ||
1099 | AliTrackPointArray *atp=tpain; | |
1100 | AliTrackPointArray *atps=NULL; | |
1101 | Int_t npts=atp->GetNPoints(); | |
1102 | ||
1103 | AliTrackPoint p; | |
1104 | // check points in specific places | |
1105 | ||
1106 | // build a new track | |
1107 | atps=new AliTrackPointArray(npts); | |
1108 | ||
1109 | Int_t npto=0; | |
1110 | for (int i=0; i<npts; i++) { | |
1111 | atp->GetPoint(p,i); | |
1112 | ||
1113 | UShort_t volid=atp->GetVolumeID()[i]; | |
1114 | Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(volid); | |
1115 | ||
b2ac0d1e | 1116 | if(index<0 || index>=2200) continue; |
eef875ec | 1117 | // dealign point |
1118 | // get MODIFIED matrix | |
1119 | TGeoHMatrix *svMatrix = m2[index]->GetSensitiveVolumeMatrix(p.GetVolumeID()); | |
1120 | //TGeoHMatrix *svOrigMatrix = mm->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID()); | |
1121 | ||
1122 | Double_t pg[3],pl[3]; | |
1123 | pg[0]=p.GetX(); | |
1124 | pg[1]=p.GetY(); | |
1125 | pg[2]=p.GetZ(); | |
1126 | if (prn) printf("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]); | |
1127 | svMatrix->MasterToLocal(pg,pl); | |
1128 | ||
1129 | // check that things went OK: local y should be 0. | |
1130 | if(TMath::Abs(pl[1])>1.e-6) { | |
1131 | printf("AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom: ERROR, local y = %f (should be zero)\n",pl[1]); | |
1132 | return -7; | |
1133 | } | |
f27a7e81 | 1134 | |
eef875ec | 1135 | if (prn) printf("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2]); |
1136 | ||
1137 | // update covariance matrix | |
1138 | TGeoHMatrix hcov; | |
1139 | Double_t hcovel[9]; | |
1140 | hcovel[0]=(Double_t)(p.GetCov()[0]); | |
1141 | hcovel[1]=(Double_t)(p.GetCov()[1]); | |
1142 | hcovel[2]=(Double_t)(p.GetCov()[2]); | |
1143 | hcovel[3]=(Double_t)(p.GetCov()[1]); | |
1144 | hcovel[4]=(Double_t)(p.GetCov()[3]); | |
1145 | hcovel[5]=(Double_t)(p.GetCov()[4]); | |
1146 | hcovel[6]=(Double_t)(p.GetCov()[2]); | |
1147 | hcovel[7]=(Double_t)(p.GetCov()[4]); | |
1148 | hcovel[8]=(Double_t)(p.GetCov()[5]); | |
1149 | hcov.SetRotation(hcovel); | |
1150 | // now rotate in local system | |
1151 | hcov.Multiply(svMatrix); | |
1152 | hcov.MultiplyLeft(&svMatrix->Inverse()); | |
1153 | // now hcov is LOCAL COVARIANCE MATRIX | |
1154 | ||
1155 | /// get original matrix of sens. vol. | |
1156 | TGeoHMatrix *svOrigMatrix = m2[index]->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID()); | |
1157 | // modify global coordinates according with pre-aligment | |
1158 | svOrigMatrix->LocalToMaster(pl,pg); | |
1159 | // now rotate in local system | |
1160 | hcov.Multiply(&svOrigMatrix->Inverse()); | |
1161 | hcov.MultiplyLeft(svOrigMatrix); | |
1162 | // hcov is back in GLOBAL RF | |
1163 | Float_t pcov[6]; | |
1164 | pcov[0]=hcov.GetRotationMatrix()[0]; | |
1165 | pcov[1]=hcov.GetRotationMatrix()[1]; | |
1166 | pcov[2]=hcov.GetRotationMatrix()[2]; | |
1167 | pcov[3]=hcov.GetRotationMatrix()[4]; | |
1168 | pcov[4]=hcov.GetRotationMatrix()[5]; | |
1169 | pcov[5]=hcov.GetRotationMatrix()[8]; | |
1170 | ||
1171 | p.SetXYZ(pg[0],pg[1],pg[2],pcov); | |
1172 | if (prn) printf("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]); | |
1173 | atps->AddPoint(npto,&p); | |
1174 | if (prn) printf("Adding point[%d] = ( %f , %f , %f ) volid = %d\n",npto,atps->GetX()[npto],atps->GetY()[npto],atps->GetZ()[npto],atps->GetVolumeID()[npto] ); | |
1175 | if (prn) p.Print(""); | |
1176 | ||
1177 | npto++; | |
1178 | } | |
1179 | ||
1180 | ||
1181 | //////////////////////////////////////////////////////////// | |
1182 | array = atps; | |
1183 | treeout->Fill(); | |
1184 | ||
1185 | delete atps; | |
1186 | atps=NULL; | |
1187 | ||
1188 | } // end loop on tracks | |
1189 | ||
1190 | pointsFile->Write(); | |
1191 | pointsFile->Close(); | |
1192 | tpainfile->Close(); | |
1193 | ||
1194 | return 0; | |
1195 | } |