]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG1/ITS/AliAlignmentDataFilterITS.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added histogram for ITS-PID cluster selection
[u/mrichter/AliRoot.git] / PWG1 / ITS / AliAlignmentDataFilterITS.cxx
CommitLineData
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
54ClassImp(AliAlignmentDataFilterITS)
55
56
b900a060 57//________________________________________________________________________
58AliAlignmentDataFilterITS::AliAlignmentDataFilterITS():
59AliAnalysisTaskSE(),
60fOnlySPDFO(kFALSE),
28647bf6 61fDownsamplelowpt(kFALSE),
b900a060 62fGeometryFileName("geometry.root"),
63fITSRecoParam(0),
64fESD(0),
65fListOfHistos(0),
66fspTree(0),
67fHistNevents(0),
68fHistNpoints(0),
69fHistPt(0),
70fHistLayer0(0),
71fHistLayer1(0),
72fHistLayer2(0),
73fHistLayer3(0),
74fHistLayer4(0),
75fHistLayer5(0),
76fntExtra(0),
77fntCosmicMatching(0)
78{
79 // Constructor
80}
81
f27a7e81 82//________________________________________________________________________
83AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
b900a060 84AliAnalysisTaskSE(name),
367c6d1f 85fOnlySPDFO(kFALSE),
28647bf6 86fDownsamplelowpt(kFALSE),
367c6d1f 87fGeometryFileName("geometry.root"),
88fITSRecoParam(0),
f27a7e81 89fESD(0),
f27a7e81 90fListOfHistos(0),
91fspTree(0),
480207c6 92fHistNevents(0),
f27a7e81 93fHistNpoints(0),
94fHistPt(0),
95fHistLayer0(0),
96fHistLayer1(0),
97fHistLayer2(0),
98fHistLayer3(0),
99fHistLayer4(0),
100fHistLayer5(0),
101fntExtra(0),
102fntCosmicMatching(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//________________________________________________________________________
115AliAlignmentDataFilterITS::~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 174void 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 259void 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());
329 cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
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//________________________________________________________________________
373void 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));
398 // Get the list of OCDB objects used for reco
399 TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
400 TIter iter2(cdbList);
401 TObjString* cdbEntry=0;
402 TString cdbEntryString;
403 while((cdbEntry =(TObjString*)(iter2.Next()))) {
404 cdbEntryString = cdbEntry->GetString();
405 if(cdbEntryString.Contains("ITS/Align/Data"))
406 itsaligndata = new TObjString(*cdbEntry);
3b9267a6 407 if(cdbEntryString.Contains("ITS/Calib/RespSDD"))
408 itscalibrespsdd = new TObjString(*cdbEntry);
ab6c74ff 409 }
410
f27a7e81 411
b900a060 412 TString triggeredClass = esd->GetFiredTriggerClasses();
367c6d1f 413 if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
414
f27a7e81 415
416 Int_t ntracks = esd->GetNumberOfTracks();
417 if(ntracks<2) return;
418
419 if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
420
421 Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
422
423 Int_t *goodtracksArray = new Int_t[ntracks];
424 Float_t *phiArray = new Float_t[ntracks];
425 Float_t *thetaArray = new Float_t[ntracks];
426 Int_t *nclsArray = new Int_t[ntracks];
427 Int_t ngt=0;
428 Int_t itrack=0;
429 for (itrack=0; itrack < ntracks; itrack++) {
430 AliESDtrack *track = esd->GetTrack(itrack);
431 if (!track) continue;
432
433
367c6d1f 434 if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
f27a7e81 435
a043746c 436 if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
437 if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
f27a7e81 438
439 Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
440 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
441
367c6d1f 442 if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
443 track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
f27a7e81 444
445 goodtracksArray[ngt] = itrack;
446 phiArray[ngt] = phi;
447 thetaArray[ngt] = theta;
448 nclsArray[ngt] = track->GetNcls(0);
449 ngt++;
450 }
451
452 if(ngt<2) {
453 delete [] goodtracksArray; goodtracksArray=0;
454 delete [] phiArray; phiArray=0;
455 delete [] thetaArray; thetaArray=0;
456 delete [] nclsArray; nclsArray=0;
457 return;
458 }
459
460 // check matching of the two tracks from the muon
461 Float_t min = 10000000.;
462 Int_t maxCls = 0;
463 Int_t good1 = -1, good2 = -1;
464 for(Int_t itr1=0; itr1<ngt-1; itr1++) {
465 for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
466 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
367c6d1f 467 if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
f27a7e81 468 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
367c6d1f 469 if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
f27a7e81 470 if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
471 maxCls = nclsArray[itr1]+nclsArray[itr2];
472 min = deltatheta+deltaphi;
473 good1 = goodtracksArray[itr1];
474 good2 = goodtracksArray[itr2];
475 } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
476 if(deltatheta+deltaphi < min) {
477 min = deltatheta+deltaphi;
478 good1 = goodtracksArray[itr1];
479 good2 = goodtracksArray[itr2];
480 }
481 }
482 }
483 }
484
485 delete [] goodtracksArray; goodtracksArray=0;
486 delete [] phiArray; phiArray=0;
487 delete [] thetaArray; thetaArray=0;
488 delete [] nclsArray; nclsArray=0;
489
490 if(good1<0) return;
491 AliDebug(2,"ok track matching");
492
493 // track1 will be the inward track (top)
494 // track2 the outward (bottom)
495 AliESDtrack *track1=0;
496 AliESDtrack *track2=0;
497 AliESDtrack *track = esd->GetTrack(good1);
498 if(track->Py()>0) {
499 track1 = esd->GetTrack(good1);
500 track2 = esd->GetTrack(good2);
501 } else {
502 track1 = esd->GetTrack(good2);
503 track2 = esd->GetTrack(good1);
504 }
505
506 AliTrackPoint point;
507 const AliTrackPointArray *array=0;
508 Int_t ipt,volId,modId,layerId,lay,lad,det;
509 Int_t jpt=0;
510 Bool_t layerOK[6][2];
511 Int_t nclsTrk[2]={0,0};
512
513 for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
514
515 for(itrack=0; itrack<2; itrack++) {
516 if(itrack==0) {
517 track = track1;
518 } else {
519 track = track2;
520 }
521 array = track->GetTrackPointArray();
522 if(!array) {
523 AliWarning("No tracks points avaialble");
524 continue;
525 }
526 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
527 array->GetPoint(point,ipt);
528 volId = point.GetVolumeID();
9f2b55aa 529 if(volId<=0) continue;
f27a7e81 530 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
9f2b55aa 531 AliDebug(2,Form("%d %d %d %f\n",ipt,layerId-1,volId,TMath::Sqrt(point.GetX()*point.GetX()+point.GetY()*point.GetY())));
f27a7e81 532 if(point.IsExtra() &&
a043746c 533 (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
9f2b55aa 534 if(layerId<1 || layerId>6) continue;
367c6d1f 535 if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
f27a7e81 536 // check minAngleWrtITSModulePlanes
537 Double_t p[3]; track->GetDirection(p);
538 TVector3 pvec(p[0],p[1],p[2]);
539 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
540 TVector3 normvec(rot[1],rot[4],rot[7]);
541 Double_t angle = pvec.Angle(normvec);
542 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
543 angle = 0.5*TMath::Pi()-angle;
367c6d1f 544 if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
f27a7e81 545 layerOK[layerId-1][itrack]=kTRUE;
546 jpt++;
547 nclsTrk[itrack]++;
548 }
549 }
550 AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
551
552 // read ITS cluster maps
553 Int_t map1[6],map2[6];
554 for(Int_t ilay=0;ilay<6;ilay++) {
555 map1[ilay]=0; map2[ilay]=0;
556 if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
557 if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
558 }
559 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()));
560 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()));
561 Int_t idx1[12],idx2[12];
562 track1->GetITSclusters(idx1);
563 track2->GetITSclusters(idx2);
564 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]));
565 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]));
566
567
367c6d1f 568 if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
f27a7e81 569 AliDebug(2,Form(" Total points %d, accepted\n",jpt));
570 fHistNpoints->Fill(jpt);
571 fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
572
573 Float_t d0z0mu[2];
574 track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
575 //printf("d0mu %f z0mu %f\n",d0z0mu[0],d0z0mu[1]);
576
44f17b29 577 vertexForTree = new AliESDVertex(*(esd->GetPrimaryVertexSPD()));
578
f27a7e81 579 Float_t dzOverlap[2];
580 Float_t curvArray[2],curverrArray[2];
581 Double_t globExtra[3],locExtra[3];
cdf8a7eb 582 if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
f27a7e81 583 arrayForTree = new AliTrackPointArray(jpt);
cdf8a7eb 584 arrayForTree->SetUniqueID(uid);
585 }
f27a7e81 586 jpt=0;
587 for(itrack=0; itrack<2; itrack++) {
588 if(itrack==0) {
589 track = track1;
590 } else {
591 track = track2;
592 }
593 curvArray[itrack] = track->GetC(esd->GetMagneticField());
594 curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
595
a043746c 596 if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
f27a7e81 597 jpt=0;
598 arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
cdf8a7eb 599 arrayForTree->SetUniqueID(uid);
f27a7e81 600 }
601 array = track->GetTrackPointArray();
602 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
603 array->GetPoint(point,ipt);
604 volId = point.GetVolumeID();
9f2b55aa 605 if(volId<=0) continue;
f27a7e81 606 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
9f2b55aa 607 if(layerId<1 || layerId>6 || !layerOK[layerId-1][itrack]) continue;
f27a7e81 608 arrayForTree->AddPoint(jpt,&point);
609 jpt++;
610 switch(layerId) {
611 case 1:
612 fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
613 break;
614 case 2:
615 fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
616 break;
617 case 3:
618 fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
619 break;
620 case 4:
621 fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
622 break;
623 case 5:
624 fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
625 break;
626 case 6:
627 fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
628 break;
629 }
b900a060 630 // Post the data for slot 2
631 if(jpt==1) PostData(2,fListOfHistos); // only if this is the first points
f27a7e81 632 if(!point.IsExtra() ||
a043746c 633 !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
f27a7e81 634 nclsTrk[itrack]--;
635 for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
636 AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
637 globExtra[0]=point.GetX();
638 globExtra[1]=point.GetY();
639 globExtra[2]=point.GetZ();
640 AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
641 //printf("%d %d %d %d %d %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
642 track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
643 AliTrackPoint pointT;
644 Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
645 for(Int_t lll=0;lll<ipt;lll++) {
646 array->GetPoint(pointT,lll);
9f2b55aa 647 if(pointT.GetVolumeID()<=0) continue;
f27a7e81 648 Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
649 if(layerIdT!=layerId) continue;
650 radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
651 radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
652 phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
653 phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
654 if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
655 thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
656 thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
657 dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
658 if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
659 dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
cdf8a7eb 660 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 661 }
662 }
663
a043746c 664 if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
f27a7e81 665 curv = curvArray[itrack];
666 curverr = curverrArray[itrack];
667 fspTree->Fill();
668 }
669 }
670
367c6d1f 671 if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
3e33351b 672 curv = 0.5*(curvArray[0]-curvArray[1]); // the "-" is because the two tracks have opposite curvature!
f27a7e81 673 curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
674 fspTree->Fill();
675 }
b900a060 676 PostData(1,fspTree);
f27a7e81 677
a043746c 678 if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return;
f27a7e81 679 // fill ntuple with track-to-track matching
680 Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;
681 Float_t d0[2],z0[2];
682 Float_t sigmad0[2],sigmaz0[2];
683 phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
684 thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
685 phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
686 thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
687 rotymu = TMath::ATan2(track1->Px(),track1->Pz());
688 rotyout = TMath::ATan2(track2->Px(),track2->Pz());
689
690 dphi = phimu - (phiout+TMath::Pi());
691 dtheta = thetamu - (TMath::Pi()-thetaout);
692 if(rotymu>0) {
693 droty = rotymu - (rotyout+TMath::Pi());
694 } else {
695 droty = rotymu - (rotyout-TMath::Pi());
696 }
697
698 Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
699
700 track1->Propagate(alpha,0.,esd->GetMagneticField());
701 track2->Propagate(alpha,0.,esd->GetMagneticField());
702 d0[0] = track1->GetY();
703 z0[0] = track1->GetZ();
704 d0[1] = track2->GetY();
705 z0[1] = track2->GetZ();
706 Float_t dxy = -(d0[0]-d0[1]);
707 Float_t dz = z0[0]-z0[1];
708 sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
709 sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
710 sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
711 sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
712 /*
713 Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
714 track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
715 track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
716 track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
717 track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
718 Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
719 Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
720 Float_t dxaty0 = x1aty0-x2aty0;
721 */
722 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]);
723
724 return;
725}
726
727//________________________________________________________________________
728void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
729{
730 // Extract ITS AliTrackPoints for Cosmics (check angular matching
731 // of top and bottom track, merge the two tracks, if requested)
732 //
733
734 // Set branch addresses for space points tree
735 AliTrackPointArray *arrayForTree=0;
44f17b29 736 AliESDVertex *vertexForTree=0;
ab6c74ff 737 Float_t curv,curverr,runNumber;
738 TObjString *itsaligndata=0;
3b9267a6 739 TObjString *itscalibrespsdd = 0;
f27a7e81 740 fspTree->SetBranchAddress("SP",&arrayForTree);
44f17b29 741 fspTree->SetBranchAddress("vertex",&vertexForTree);
f27a7e81 742 fspTree->SetBranchAddress("curv",&curv);
743 fspTree->SetBranchAddress("curverr",&curverr);
ab6c74ff 744 fspTree->SetBranchAddress("run",&runNumber);
745 fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
3b9267a6 746 fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
ab6c74ff 747
748
749 runNumber = (Float_t)esd->GetRunNumber();
b900a060 750 Int_t uid=20000+esd->GetEventNumberInFile();
ab6c74ff 751
752 TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
753 // Get the list of OCDB objects used for reco
754 TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
755 TIter iter2(cdbList);
756 TObjString* cdbEntry=0;
757 TString cdbEntryString;
758 while((cdbEntry =(TObjString*)(iter2.Next()))) {
759 cdbEntryString = cdbEntry->GetString();
760 if(cdbEntryString.Contains("ITS/Align/Data"))
761 itsaligndata = new TObjString(*cdbEntry);
3b9267a6 762 if(cdbEntryString.Contains("ITS/Calib/RespSDD"))
763 itscalibrespsdd = new TObjString(*cdbEntry);
ab6c74ff 764 }
f27a7e81 765
766 Int_t ntracks = esd->GetNumberOfTracks();
767
768 if(ntracks==0) return;
769
44f17b29 770 const AliESDVertex *vertexTracks = esd->GetPrimaryVertexTracks();
771 if(!vertexTracks);
772 if(vertexTracks->GetNContributors()<=0) return;
f27a7e81 773
44f17b29 774 Double_t vtxpos[3]; vertexTracks->GetXYZ(vtxpos);
f27a7e81 775
776 Int_t ncls=0;
777 Double_t pt=-10000.;
778 Double_t d0z0[2],covd0z0[3];
779 const AliTrackPointArray *array = 0;
780
781 for (Int_t itrack=0; itrack < ntracks; itrack++) {
782 AliESDtrack * track = esd->GetTrack(itrack);
783 if (!track) continue;
784
28647bf6 785 if(fDownsamplelowpt && TMath::Abs(esd->GetMagneticField())>0.01 &&
786 track->Pt()<gRandom->Rndm()) continue;
787
367c6d1f 788 if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
f27a7e81 789
a043746c 790 if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
791 if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
f27a7e81 792
367c6d1f 793 if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
794 track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
f27a7e81 795
796 pt = track->Pt();
797 ncls = track->GetNcls(0);
798 Double_t maxd=10000.;
44f17b29 799 track->PropagateToDCA(vertexTracks,esd->GetMagneticField(),maxd,d0z0,covd0z0);
f27a7e81 800
801 // read ITS cluster map
802 Int_t map[6];
803 for(Int_t ilay=0;ilay<6;ilay++) {
804 map[ilay]=0;
805 if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
806 }
807 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()));
808 Int_t idx[12];
809 track->GetITSclusters(idx);
810 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]));
811
812
813 AliTrackPoint point;
814 Int_t ipt,volId,modId,layerId,lay,lad,det;
815 Int_t jpt=0;
816 Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
817
818 array = track->GetTrackPointArray();
28647bf6 819 if(!array) {printf("no track points\n"); continue;}
f27a7e81 820 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
821 array->GetPoint(point,ipt);
822 volId = point.GetVolumeID();
9f2b55aa 823 if(volId<=0) continue;
f27a7e81 824 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
825 if(layerId<1 || layerId>6) continue;
826 if(point.IsExtra() &&
a043746c 827 (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
f27a7e81 828 layerOK[layerId-1]=kTRUE;
829 jpt++;
830 }
831
367c6d1f 832 if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
f27a7e81 833
834 fHistNpoints->Fill(jpt);
835 fHistPt->Fill(pt);
b900a060 836 PostData(2,fListOfHistos);
f27a7e81 837
838 Float_t dzOverlap[2];
839 Double_t globExtra[3],locExtra[3];
840 arrayForTree = new AliTrackPointArray(jpt);
cdf8a7eb 841 arrayForTree->SetUniqueID(uid);
f27a7e81 842 jpt=0;
843 array = track->GetTrackPointArray();
844 if(!array) continue;
845 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
846 array->GetPoint(point,ipt);
847 volId = point.GetVolumeID();
848 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
849 if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
44c3f5f6 850 arrayForTree->AddPoint(jpt,&point);
851 switch(layerId) {
852 case 1:
853 fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
854 break;
855 case 2:
856 fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
857 break;
858 case 3:
859 fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
860 break;
861 case 4:
862 fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
863 break;
864 case 5:
865 fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
866 break;
867 case 6:
868 fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
869 break;
870 }
871 jpt++;
f27a7e81 872 if(!point.IsExtra() ||
a043746c 873 !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
f27a7e81 874 ncls--;
875 for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
876 AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
877 globExtra[0]=point.GetX();
878 globExtra[1]=point.GetY();
879 globExtra[2]=point.GetZ();
880 AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
881 track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
882 AliTrackPoint pointT;
883 Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
884 for(Int_t lll=0;lll<ipt;lll++) {
885 array->GetPoint(pointT,lll);
886 Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
887 if(layerIdT!=layerId) continue;
888 radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
889 radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
890 phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
891 phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
892 if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
893 thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
894 thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
895 dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
896 if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
897 dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
cdf8a7eb 898 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 899 }
f27a7e81 900 }
901
902 curv = track->GetC(esd->GetMagneticField());
903 curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
904
44f17b29 905 vertexForTree = new AliESDVertex(*vertexTracks);
906 if(vertexTracks->UsesTrack(track->GetID())) {
907 vertexForTree->SetID(1);
908 } else {
909 vertexForTree->SetID(0);
910 }
911
f27a7e81 912 fspTree->Fill();
913
914 } // end of tracks loop
915
b900a060 916 PostData(1,fspTree);
f27a7e81 917
918 return;
919}
920
921//________________________________________________________________________
922void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
923{
924 // Terminate analysis
925 //
926 AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
927
81e84469 928 fspTree = dynamic_cast<TTree*> (GetOutputData(1));
f27a7e81 929 if (!fspTree) {
930 printf("ERROR: fspTree not available\n");
931 return;
932 }
933
81e84469 934 fListOfHistos = dynamic_cast<TList*> (GetOutputData(2));
f27a7e81 935 if (!fListOfHistos) {
936 printf("ERROR: fListOfHistos not available\n");
937 return;
938 }
939
a043746c 940 fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents"));
f27a7e81 941 fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
942 fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
943 fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
944 fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
945 fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
946 fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
947 fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
948 fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
949 fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
950 fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
951
952
953
954 return;
955}
eef875ec 956//-------------------------------------------------------------------------------
957const AliITSRecoParam *AliAlignmentDataFilterITS::GetRecoParam() const
958{
959 //
960 // Return the ITSRecoParam object
961 //
962 if(AliITSReconstructor::GetRecoParam()) {
963 return AliITSReconstructor::GetRecoParam();
964 } else if(fITSRecoParam) {
965 return fITSRecoParam;
966 } else return NULL;
967}
968//--------------------------------------------------------------------------------
969Int_t AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom(Char_t *fin,
970 Char_t *fout,
971 Char_t *fmis,
972 Char_t *fgeo,
973 Bool_t prn)
974{
975 //
976 // Convert AliTrackPoints in fin, reconstructed with fmis, back
977 // to ideal geometry
978 //
979 // M. Lunardon
980 //
981
982
983 TGeoHMatrix deltahm;
eef875ec 984
985 // Load geometry
986 if (gSystem->AccessPathName(fgeo)) {
987 printf("couldn't find geometry file %s - skipping...\n",fmis);
988 return -1;
989 }
990
991 TFile *geofile=TFile::Open(fgeo);
992 TGeoManager *fgGeometry=NULL;
993
994 fgGeometry=(TGeoManager*)geofile->Get("ALICE");
995
996 if (!fgGeometry)
997 fgGeometry=(TGeoManager*)geofile->Get("Geometry");
998
999 if (!fgGeometry) {
1000 AliCDBEntry *entry = (AliCDBEntry*)geofile->Get("AliCDBEntry");
1001 if (entry)
1002 fgGeometry = (TGeoManager*)entry->GetObject();
1003 }
1004
1005 if (!fgGeometry) return -1;
1006 AliGeomManager::SetGeometry(fgGeometry);
1007 if(!AliGeomManager::GetGeometry()) return -1;
1008
1009
1010 // open alignment file
1011 if (gSystem->AccessPathName(fmis)) {
1012 printf("couldn't open alignment file %s - skipping...\n",fmis);
1013 return -2;
1014 }
1015 TFile *pref = TFile::Open(fmis);
1016 if (!pref->IsOpen()) return -2;
1017
1018
1019 /// apply alignment to ideal geometry
1020 TClonesArray *prea=(TClonesArray*)pref->Get("ITSAlignObjs");
1021 if (!prea) {
1022 if (pref->Get("AliCDBEntry"))
1023 prea = (TClonesArray*) ((AliCDBEntry*)pref->Get("AliCDBEntry"))->GetObject();
1024 }
1025 if (!prea) return -3;
1026 Int_t nprea=prea->GetEntriesFast();
1027 printf("Array of input misalignments with %d entries\n",nprea);
1028 AliGeomManager::ApplyAlignObjsToGeom(*prea); // apply all levels of objs
1029
1030 AliTrackPointArray *tpain=NULL;
1031 TFile *tpainfile=NULL;
1032 TTree *treein=NULL;
1033 AliTrackPoint point;
1034 AliITSAlignMille2Module *m2[2200];
1035 for (Int_t i=0; i<2198; i++)
1036 m2[i]=new AliITSAlignMille2Module(AliITSAlignMille2Module::GetVolumeIDFromIndex(i));
1037
1038 // open input file
1039 if (gSystem->AccessPathName(fin)) {
1040 printf("couldn't open file %s - skipping...\n",fin);
1041 return -4;
1042 }
1043 tpainfile = TFile::Open(fin);
1044 if (!tpainfile->IsOpen()) return -4;
1045
1046 treein=(TTree*)tpainfile->Get("spTree");
1047 if (!treein) return -5;
1048 Float_t curv,curverr,runNumber;
1049 TObjString *itsaligndata=0;
1050 TObjString *itscalibrespsdd = 0;
1051 treein->SetBranchAddress("SP", &tpain);
1052 treein->SetBranchAddress("curv", &curv);
1053 treein->SetBranchAddress("curverr", &curverr);
1054 treein->SetBranchAddress("run",&runNumber);
1055 treein->SetBranchAddress("ITSAlignData",&itsaligndata);
1056 treein->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
1057
1058 int ntrks=treein->GetEntries();
1059 printf("Reading %d tracks from %s\n",ntrks,fin);
1060
1061
1062 // open output file
1063 TFile *pointsFile = TFile::Open(fout,"RECREATE");
1064 if (!pointsFile || !pointsFile->IsOpen()) {
1065 printf("Can't open output file %s !",fout);
1066 return -6;
1067 }
1068 AliTrackPointArray *array = new AliTrackPointArray();
1069
1070 // new!
1071 TTree *treeout=(TTree*)treein->Clone("spTree");
1072 treeout->Reset();
1073 treeout->SetBranchAddress("SP", &array);
1074 treeout->SetBranchAddress("curv", &curv);
1075 treeout->SetBranchAddress("curverr", &curverr);
1076 treeout->SetBranchAddress("run",&runNumber);
1077 treeout->SetBranchAddress("ITSAlignData",&itsaligndata);
1078 treeout->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
1079
1080 // tracks main loop
1081 for (Int_t it=0; it<ntrks; it++) {
1082 if (!(it%5000) ) printf("...processing track n. %d\n",it);
1083
1084 treein->GetEvent(it);
1085
1086 //////////////////////////////
1087
1088 AliTrackPointArray *atp=tpain;
1089 AliTrackPointArray *atps=NULL;
1090 Int_t npts=atp->GetNPoints();
1091
1092 AliTrackPoint p;
1093 // check points in specific places
1094
1095 // build a new track
1096 atps=new AliTrackPointArray(npts);
1097
1098 Int_t npto=0;
1099 for (int i=0; i<npts; i++) {
1100 atp->GetPoint(p,i);
1101
1102 UShort_t volid=atp->GetVolumeID()[i];
1103 Int_t index=AliITSAlignMille2Module::GetIndexFromVolumeID(volid);
1104
1105
1106 // dealign point
1107 // get MODIFIED matrix
1108 TGeoHMatrix *svMatrix = m2[index]->GetSensitiveVolumeMatrix(p.GetVolumeID());
1109 //TGeoHMatrix *svOrigMatrix = mm->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1110
1111 Double_t pg[3],pl[3];
1112 pg[0]=p.GetX();
1113 pg[1]=p.GetY();
1114 pg[2]=p.GetZ();
1115 if (prn) printf("Global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]);
1116 svMatrix->MasterToLocal(pg,pl);
1117
1118 // check that things went OK: local y should be 0.
1119 if(TMath::Abs(pl[1])>1.e-6) {
1120 printf("AliAlignmentDataFilterITS::WriteTrackPointsInIdealGeom: ERROR, local y = %f (should be zero)\n",pl[1]);
1121 return -7;
1122 }
f27a7e81 1123
eef875ec 1124 if (prn) printf("Local coordinates of measured point : X=%f Y=%f Z=%f \n",pl[0],pl[1],pl[2]);
1125
1126 // update covariance matrix
1127 TGeoHMatrix hcov;
1128 Double_t hcovel[9];
1129 hcovel[0]=(Double_t)(p.GetCov()[0]);
1130 hcovel[1]=(Double_t)(p.GetCov()[1]);
1131 hcovel[2]=(Double_t)(p.GetCov()[2]);
1132 hcovel[3]=(Double_t)(p.GetCov()[1]);
1133 hcovel[4]=(Double_t)(p.GetCov()[3]);
1134 hcovel[5]=(Double_t)(p.GetCov()[4]);
1135 hcovel[6]=(Double_t)(p.GetCov()[2]);
1136 hcovel[7]=(Double_t)(p.GetCov()[4]);
1137 hcovel[8]=(Double_t)(p.GetCov()[5]);
1138 hcov.SetRotation(hcovel);
1139 // now rotate in local system
1140 hcov.Multiply(svMatrix);
1141 hcov.MultiplyLeft(&svMatrix->Inverse());
1142 // now hcov is LOCAL COVARIANCE MATRIX
1143
1144 /// get original matrix of sens. vol.
1145 TGeoHMatrix *svOrigMatrix = m2[index]->GetSensitiveVolumeOrigGlobalMatrix(p.GetVolumeID());
1146 // modify global coordinates according with pre-aligment
1147 svOrigMatrix->LocalToMaster(pl,pg);
1148 // now rotate in local system
1149 hcov.Multiply(&svOrigMatrix->Inverse());
1150 hcov.MultiplyLeft(svOrigMatrix);
1151 // hcov is back in GLOBAL RF
1152 Float_t pcov[6];
1153 pcov[0]=hcov.GetRotationMatrix()[0];
1154 pcov[1]=hcov.GetRotationMatrix()[1];
1155 pcov[2]=hcov.GetRotationMatrix()[2];
1156 pcov[3]=hcov.GetRotationMatrix()[4];
1157 pcov[4]=hcov.GetRotationMatrix()[5];
1158 pcov[5]=hcov.GetRotationMatrix()[8];
1159
1160 p.SetXYZ(pg[0],pg[1],pg[2],pcov);
1161 if (prn) printf("New global coordinates of measured point : X=%f Y=%f Z=%f \n",pg[0],pg[1],pg[2]);
1162 atps->AddPoint(npto,&p);
1163 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] );
1164 if (prn) p.Print("");
1165
1166 npto++;
1167 }
1168
1169
1170 ////////////////////////////////////////////////////////////
1171 array = atps;
1172 treeout->Fill();
1173
1174 delete atps;
1175 atps=NULL;
1176
1177 } // end loop on tracks
1178
1179 pointsFile->Write();
1180 pointsFile->Close();
1181 tpainfile->Close();
1182
1183 return 0;
1184}