]> 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
Fix
[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>
25#include <TChain.h>
26#include <TNtuple.h>
27#include <TList.h>
28#include <TH1F.h>
29#include <TH2F.h>
cf53e6db 30#include <TMap.h>
f27a7e81 31#include <TVector3.h>
32#include <TGeoManager.h>
33
34#include "AliLog.h"
35#include "AliGeomManager.h"
a043746c 36#include "AliITSReconstructor.h"
f27a7e81 37#include "AliITSgeomTGeo.h"
38#include "AliTrackPointArray.h"
39#include "AliESDInputHandler.h"
40#include "AliESDVertex.h"
41#include "AliESDtrack.h"
42#include "AliESDEvent.h"
43#include "AliESDfriend.h"
44#include "AliAnalysisTask.h"
45#include "AliAnalysisManager.h"
46#include "AliAlignmentDataFilterITS.h"
47
48
49ClassImp(AliAlignmentDataFilterITS)
50
51
52//________________________________________________________________________
53AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
54AliAnalysisTask(name,"task"),
367c6d1f 55fOnlySPDFO(kFALSE),
56fGeometryFileName("geometry.root"),
57fITSRecoParam(0),
f27a7e81 58fESD(0),
59fESDfriend(0),
60fListOfHistos(0),
61fspTree(0),
480207c6 62fHistNevents(0),
f27a7e81 63fHistNpoints(0),
64fHistPt(0),
65fHistLayer0(0),
66fHistLayer1(0),
67fHistLayer2(0),
68fHistLayer3(0),
69fHistLayer4(0),
70fHistLayer5(0),
71fntExtra(0),
72fntCosmicMatching(0)
73{
74 // Constructor
75
76 // Define input and output slots here
77 // Input slot #0 works with a TChain
78 DefineInput(0, TChain::Class());
79 // Output slot #0 writes into a TTree
80 DefineOutput(0,TTree::Class()); //My private output
81 // Output slot #1 writes into a TList
82 DefineOutput(1,TList::Class()); //My private output
83}
84
85//________________________________________________________________________
86AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS()
87{
88 // Destructor
89 if (fListOfHistos) {
90 delete fListOfHistos;
91 fListOfHistos = 0;
92 }
93 if (fspTree) {
94 delete fspTree;
95 fspTree = 0;
96 }
480207c6 97 if (fHistNevents) {
98 delete fHistNevents;
99 fHistNevents = 0;
100 }
a043746c 101 if (fHistNpoints) {
102 delete fHistNpoints;
103 fHistNpoints = 0;
104 }
f27a7e81 105 if (fHistPt) {
106 delete fHistPt;
107 fHistPt = 0;
108 }
109 if (fHistLayer0) {
110 delete fHistLayer0;
111 fHistLayer0 = 0;
112 }
113 if (fHistLayer1) {
114 delete fHistLayer1;
115 fHistLayer1 = 0;
116 }
117 if (fHistLayer2) {
118 delete fHistLayer2;
119 fHistLayer2 = 0;
120 }
121 if (fHistLayer3) {
122 delete fHistLayer3;
123 fHistLayer3 = 0;
124 }
125 if (fHistLayer4) {
126 delete fHistLayer4;
127 fHistLayer4 = 0;
128 }
129 if (fHistLayer5) {
130 delete fHistLayer5;
131 fHistLayer5 = 0;
132 }
133 if (fntExtra) {
134 delete fntExtra;
135 fntExtra = 0;
136 }
137 if (fntCosmicMatching) {
138 delete fntCosmicMatching;
139 fntCosmicMatching = 0;
140 }
141}
a043746c 142
f27a7e81 143//________________________________________________________________________
144void AliAlignmentDataFilterITS::ConnectInputData(Option_t *)
145{
146 // Connect ESD
147 // Called once
148
149 TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
150 if(!tree) {
151 printf("ERROR: Could not read chain from input slot 0\n");
152 } else {
cf53e6db 153 // Get the OCDB path and the list of OCDB objects used for reco
154 TMap *cdbMap = (TMap*)(tree->GetTree()->GetUserInfo())->FindObject("cdbMap");
155 TList *cdbList = (TList*)(tree->GetTree()->GetUserInfo())->FindObject("cdbList");
156
157 //cdbList->Print();
158 // write the list to the user info of the output tree
159 if(!fspTree) {
160 printf("ERROR: fspTree does not exist\n");
161 } else {
162 TMap *cdbMapCopy = new TMap(cdbMap->GetEntries());
163 cdbMapCopy->SetOwner(1);
164 cdbMapCopy->SetName("cdbMap");
165 TIter iter1(cdbMap->GetTable());
166
167 TPair* pair = 0;
168 while((pair = dynamic_cast<TPair*> (iter1.Next()))){
169 TObjString* keyStr = dynamic_cast<TObjString*> (pair->Key());
170 TObjString* valStr = dynamic_cast<TObjString*> (pair->Value());
171 cdbMapCopy->Add(new TObjString(keyStr->GetName()), new TObjString(valStr->GetName()));
172 }
173
174 TList *cdbListCopy = new TList();
175 cdbListCopy->SetOwner(1);
176 cdbListCopy->SetName("cdbList");
177
178 TIter iter2(cdbList);
179
180 TObjString* cdbEntry=0;
181 while((cdbEntry =(TObjString*)(iter2.Next()))) {
182 cdbListCopy->Add(new TObjString(*cdbEntry));
183 }
184 cdbListCopy->Print();
185
186
187 fspTree->GetUserInfo()->Add(cdbMapCopy);
188 fspTree->GetUserInfo()->Add(cdbListCopy);
189 }
f27a7e81 190
cf53e6db 191 // Disable all branches and enable only the needed ones
f27a7e81 192 tree->SetBranchStatus("fTriggerMask", 1);
193 tree->SetBranchStatus("fSPDVertex*", 1);
f27a7e81 194 tree->SetBranchStatus("ESDfriend*", 1);
195 tree->SetBranchAddress("ESDfriend.",&fESDfriend);
f27a7e81 196 tree->SetBranchStatus("fSPDMult*", 1);
197
198 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
199
200 if(!esdH) {
201 printf("ERROR: Could not get ESDInputHandler\n");
202 } else {
203 fESD = esdH->GetEvent();
204 }
205 }
206
207 return;
208}
209
210//________________________________________________________________________
211void AliAlignmentDataFilterITS::Init()
212{
213 // Initialization
214
215 return;
216}
217
218//________________________________________________________________________
219void AliAlignmentDataFilterITS::CreateOutputObjects()
220{
221 // Create the output container
222 //
c1605e49 223
224 // load the geometry
225 if(!gGeoManager) {
226 printf("AliAlignmentDataFilterITS::CreateOutputObjects(): loading geometry from %s\n",fGeometryFileName.Data());
227 AliGeomManager::LoadGeometry(fGeometryFileName.Data());
228 if(!gGeoManager) {
229 printf("AliAlignmentDataFilterITS::CreateOutputObjects(): no geometry loaded \n");
230 return;
231 }
232 }
f27a7e81 233
234 // Several histograms are more conveniently managed in a TList
235 fListOfHistos = new TList();
236 fListOfHistos->SetOwner();
a043746c 237
480207c6 238 fHistNevents = new TH1F("fHistNevents", "Number of processed events; N events; bin",5,-0.5,4.5);
239 fHistNevents->Sumw2();
240 fHistNevents->SetMinimum(0);
241 fListOfHistos->Add(fHistNevents);
f27a7e81 242
243 fHistNpoints = new TH1F("fHistNpoints", "Number of AliTrackPoints per track; N points; tracks",25,-0.5,24.5);
244 fHistNpoints->Sumw2();
245 fHistNpoints->SetMinimum(0);
246 fListOfHistos->Add(fHistNpoints);
247
248 fHistPt = new TH1F("fHistPt", "p_{t} of tracks; p_{t} [GeV/c]; tracks",100,0,50);
249 fHistPt->Sumw2();
250 fHistPt->SetMinimum(0);
251 fListOfHistos->Add(fHistPt);
252
253
254 Float_t zmax=14.;
255 Int_t nbinsphi=20,nbinsz=4;
256 fHistLayer0 = new TH2F("fHistLayer0","Points in layer inner SPD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
257 fListOfHistos->Add(fHistLayer0);
258 zmax=14.;
259 nbinsphi=40;nbinsz=4;
260 fHistLayer1 = new TH2F("fHistLayer1","Points in layer outer SPD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
261 fListOfHistos->Add(fHistLayer1);
262 zmax=22.;
263 nbinsphi=14;nbinsz=6;
264 fHistLayer2 = new TH2F("fHistLayer2","Points in layer inner SDD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
265 fListOfHistos->Add(fHistLayer2);
266 zmax=29.5;
267 nbinsphi=22;nbinsz=8;
268 fHistLayer3 = new TH2F("fHistLayer3","Points in layer outer SDD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
269 fListOfHistos->Add(fHistLayer3);
270 zmax=45.;
271 nbinsphi=34;nbinsz=23;
272 fHistLayer4 = new TH2F("fHistLayer4","Points in layer inner SSD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
273 fListOfHistos->Add(fHistLayer4);
274 zmax=51.;
275 nbinsphi=38;nbinsz=26;
276 fHistLayer5 = new TH2F("fHistLayer5","Points in layer outer SSD; global #phi; global z [cm]",nbinsphi,-3.14,3.14,nbinsz,-zmax,zmax);
277 fListOfHistos->Add(fHistLayer5);
278
279
280 fntExtra = new TNtuple("fntExtra","extra clusters in ITS","ncls:layer:ladder:volid:phi:x:y:z:xloc:zloc:dxy:dz:d0mu:z0mu");
281 fListOfHistos->Add(fntExtra);
282
283 fntCosmicMatching = new TNtuple("fntCosmicMatching","cosmic tracks matching in ITS","ncls1:ncls2:pt1:pt2:sigmad01:sigmad02:sigmaz01:sigmaz02:dxy:dz:phimu:thetamu:d0mu:z0mu");
284 fListOfHistos->Add(fntCosmicMatching);
285
3acc122a 286 fspTree = new TTree("spTree","Tree with ITS track points");
f27a7e81 287 const AliTrackPointArray *array = 0;
ab6c74ff 288 Float_t curv,curverr,runNumber;
289 const TObjString *itsaligndata = 0;
3b9267a6 290 const TObjString *itscalibrespsdd = 0;
f27a7e81 291 fspTree->Branch("SP","AliTrackPointArray",&array);
292 fspTree->Branch("curv",&curv);
293 fspTree->Branch("curverr",&curverr);
ab6c74ff 294 fspTree->Branch("run",&runNumber);
295 fspTree->Branch("ITSAlignData",&itsaligndata);
3b9267a6 296 fspTree->Branch("ITSCalibRespSDD",&itscalibrespsdd);
ab6c74ff 297
f27a7e81 298
299 return;
300}
301
302//________________________________________________________________________
303void AliAlignmentDataFilterITS::Exec(Option_t */*option*/)
304{
305 // Execute analysis for current event:
306 // write ITS AliTrackPoints for selected tracks to fspTree
a043746c 307
c1605e49 308 // check the geometry
309 if(!gGeoManager) {
310 printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
311 return;
367c6d1f 312 }
313
314 // check if we have AliITSRecoParam
a043746c 315 if(!GetRecoParam()) {
367c6d1f 316 if(!fITSRecoParam) {
317 printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n");
318 return;
319 }
f27a7e81 320 }
321
a043746c 322
f27a7e81 323 if(!fESD) {
324 printf("AliAlignmentDataFilterITS::Exec(): no ESD \n");
325 return;
326 }
327 if(!fESDfriend) {
328 printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n");
329 return;
330 }
331 // attach ESDfriend
332 fESD->SetESDfriend(fESDfriend);
333
480207c6 334 // Post the data for slot 0
335 fHistNevents->Fill(0);
336
cf53e6db 337
338 // write field value to spTree UserInfo
339 if(!((fspTree->GetUserInfo())->FindObject("BzkGauss"))) {
340 Double_t bz=fESD->GetMagneticField();
341 TString bzString; bzString+=bz;
342 TObjString *bzObjString = new TObjString(bzString);
343 TList *bzList = new TList();
344 bzList->SetOwner(1);
345 bzList->SetName("BzkGauss");
346 bzList->Add(bzObjString);
347 fspTree->GetUserInfo()->Add(bzList);
348 }
349
350
f27a7e81 351 // Process event as Cosmic or Collision
352 //if(esd->GetEventType()== ???? ) {
353 printf("AliAlignmentDataFilterITS::Exec(): MOVE ASAP TO esd->GetEventType() !\n");
367c6d1f 354 if(GetRecoParam()->GetAlignFilterCosmics()) {
f27a7e81 355 FilterCosmic(fESD);
356 } else {
357 FilterCollision(fESD);
358 }
359
480207c6 360 PostData(1,fListOfHistos);
361
f27a7e81 362 return;
363}
364
365//________________________________________________________________________
366void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd)
367{
368 // Extract ITS AliTrackPoints for Cosmics (check angular matching
369 // of top and bottom track, merge the two tracks, if requested)
370 //
371
372 // Set branch addresses for space points tree
373 AliTrackPointArray *arrayForTree=0;
ab6c74ff 374 Float_t curv,curverr,runNumber;
375 TObjString *itsaligndata=0;
3b9267a6 376 TObjString *itscalibrespsdd = 0;
f27a7e81 377 fspTree->SetBranchAddress("SP",&arrayForTree);
378 fspTree->SetBranchAddress("curv",&curv);
379 fspTree->SetBranchAddress("curverr",&curverr);
ab6c74ff 380 fspTree->SetBranchAddress("run",&runNumber);
381 fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
3b9267a6 382 fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
ab6c74ff 383
384
385 runNumber = (Float_t)esd->GetRunNumber();
386
387 TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
388 // Get the list of OCDB objects used for reco
389 TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
390 TIter iter2(cdbList);
391 TObjString* cdbEntry=0;
392 TString cdbEntryString;
393 while((cdbEntry =(TObjString*)(iter2.Next()))) {
394 cdbEntryString = cdbEntry->GetString();
395 if(cdbEntryString.Contains("ITS/Align/Data"))
396 itsaligndata = new TObjString(*cdbEntry);
3b9267a6 397 if(cdbEntryString.Contains("ITS/Calib/RespSDD"))
398 itscalibrespsdd = new TObjString(*cdbEntry);
ab6c74ff 399 }
400
f27a7e81 401
367c6d1f 402 TString triggeredClass = fESD->GetFiredTriggerClasses();
403 if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
404
f27a7e81 405
406 Int_t ntracks = esd->GetNumberOfTracks();
407 if(ntracks<2) return;
408
409 if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
410
411 Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
412
413 Int_t *goodtracksArray = new Int_t[ntracks];
414 Float_t *phiArray = new Float_t[ntracks];
415 Float_t *thetaArray = new Float_t[ntracks];
416 Int_t *nclsArray = new Int_t[ntracks];
417 Int_t ngt=0;
418 Int_t itrack=0;
419 for (itrack=0; itrack < ntracks; itrack++) {
420 AliESDtrack *track = esd->GetTrack(itrack);
421 if (!track) continue;
422
423
367c6d1f 424 if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
f27a7e81 425
a043746c 426 if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
427 if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
f27a7e81 428
429 Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
430 Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
431
367c6d1f 432 if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
433 track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
f27a7e81 434
435 goodtracksArray[ngt] = itrack;
436 phiArray[ngt] = phi;
437 thetaArray[ngt] = theta;
438 nclsArray[ngt] = track->GetNcls(0);
439 ngt++;
440 }
441
442 if(ngt<2) {
443 delete [] goodtracksArray; goodtracksArray=0;
444 delete [] phiArray; phiArray=0;
445 delete [] thetaArray; thetaArray=0;
446 delete [] nclsArray; nclsArray=0;
447 return;
448 }
449
450 // check matching of the two tracks from the muon
451 Float_t min = 10000000.;
452 Int_t maxCls = 0;
453 Int_t good1 = -1, good2 = -1;
454 for(Int_t itr1=0; itr1<ngt-1; itr1++) {
455 for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
456 Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
367c6d1f 457 if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
f27a7e81 458 Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
367c6d1f 459 if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
f27a7e81 460 if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
461 maxCls = nclsArray[itr1]+nclsArray[itr2];
462 min = deltatheta+deltaphi;
463 good1 = goodtracksArray[itr1];
464 good2 = goodtracksArray[itr2];
465 } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
466 if(deltatheta+deltaphi < min) {
467 min = deltatheta+deltaphi;
468 good1 = goodtracksArray[itr1];
469 good2 = goodtracksArray[itr2];
470 }
471 }
472 }
473 }
474
475 delete [] goodtracksArray; goodtracksArray=0;
476 delete [] phiArray; phiArray=0;
477 delete [] thetaArray; thetaArray=0;
478 delete [] nclsArray; nclsArray=0;
479
480 if(good1<0) return;
481 AliDebug(2,"ok track matching");
482
483 // track1 will be the inward track (top)
484 // track2 the outward (bottom)
485 AliESDtrack *track1=0;
486 AliESDtrack *track2=0;
487 AliESDtrack *track = esd->GetTrack(good1);
488 if(track->Py()>0) {
489 track1 = esd->GetTrack(good1);
490 track2 = esd->GetTrack(good2);
491 } else {
492 track1 = esd->GetTrack(good2);
493 track2 = esd->GetTrack(good1);
494 }
495
496 AliTrackPoint point;
497 const AliTrackPointArray *array=0;
498 Int_t ipt,volId,modId,layerId,lay,lad,det;
499 Int_t jpt=0;
500 Bool_t layerOK[6][2];
501 Int_t nclsTrk[2]={0,0};
502
503 for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
504
505 for(itrack=0; itrack<2; itrack++) {
506 if(itrack==0) {
507 track = track1;
508 } else {
509 track = track2;
510 }
511 array = track->GetTrackPointArray();
512 if(!array) {
513 AliWarning("No tracks points avaialble");
514 continue;
515 }
516 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
517 array->GetPoint(point,ipt);
518 volId = point.GetVolumeID();
519 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
520 AliDebug(2,Form("%d %d\n",ipt,layerId-1));
521 if(point.IsExtra() &&
a043746c 522 (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
f27a7e81 523 if(layerId>6) continue;
367c6d1f 524 if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
f27a7e81 525 // check minAngleWrtITSModulePlanes
526 Double_t p[3]; track->GetDirection(p);
527 TVector3 pvec(p[0],p[1],p[2]);
528 Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
529 TVector3 normvec(rot[1],rot[4],rot[7]);
530 Double_t angle = pvec.Angle(normvec);
531 if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
532 angle = 0.5*TMath::Pi()-angle;
367c6d1f 533 if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
f27a7e81 534 layerOK[layerId-1][itrack]=kTRUE;
535 jpt++;
536 nclsTrk[itrack]++;
537 }
538 }
539 AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
540
541 // read ITS cluster maps
542 Int_t map1[6],map2[6];
543 for(Int_t ilay=0;ilay<6;ilay++) {
544 map1[ilay]=0; map2[ilay]=0;
545 if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
546 if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
547 }
548 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()));
549 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()));
550 Int_t idx1[12],idx2[12];
551 track1->GetITSclusters(idx1);
552 track2->GetITSclusters(idx2);
553 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]));
554 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]));
555
556
367c6d1f 557 if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
f27a7e81 558 AliDebug(2,Form(" Total points %d, accepted\n",jpt));
559 fHistNpoints->Fill(jpt);
560 fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
561
562 Float_t d0z0mu[2];
563 track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
564 //printf("d0mu %f z0mu %f\n",d0z0mu[0],d0z0mu[1]);
565
566 Float_t dzOverlap[2];
567 Float_t curvArray[2],curverrArray[2];
568 Double_t globExtra[3],locExtra[3];
367c6d1f 569 if(GetRecoParam()->GetAlignFilterCosmicMergeTracks())
f27a7e81 570 arrayForTree = new AliTrackPointArray(jpt);
571
572 jpt=0;
573 for(itrack=0; itrack<2; itrack++) {
574 if(itrack==0) {
575 track = track1;
576 } else {
577 track = track2;
578 }
579 curvArray[itrack] = track->GetC(esd->GetMagneticField());
580 curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
581
a043746c 582 if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
f27a7e81 583 jpt=0;
584 arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
585 }
586 array = track->GetTrackPointArray();
587 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
588 array->GetPoint(point,ipt);
589 volId = point.GetVolumeID();
590 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
591 if(layerId>6 || !layerOK[layerId-1][itrack]) continue;
592 arrayForTree->AddPoint(jpt,&point);
593 jpt++;
594 switch(layerId) {
595 case 1:
596 fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
597 break;
598 case 2:
599 fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
600 break;
601 case 3:
602 fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
603 break;
604 case 4:
605 fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
606 break;
607 case 5:
608 fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
609 break;
610 case 6:
611 fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
612 break;
613 }
614 // Post the data for slot 0
615 if(jpt==1) PostData(1,fListOfHistos); // only if this is the first points
616 if(!point.IsExtra() ||
a043746c 617 !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
f27a7e81 618 nclsTrk[itrack]--;
619 for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
620 AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
621 globExtra[0]=point.GetX();
622 globExtra[1]=point.GetY();
623 globExtra[2]=point.GetZ();
624 AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
625 //printf("%d %d %d %d %d %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
626 track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
627 AliTrackPoint pointT;
628 Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
629 for(Int_t lll=0;lll<ipt;lll++) {
630 array->GetPoint(pointT,lll);
631 Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
632 if(layerIdT!=layerId) continue;
633 radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
634 radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
635 phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
636 phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
637 if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
638 thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
639 thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
640 dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
641 if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
642 dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
643 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]);
644 }
645 }
646
a043746c 647 if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
f27a7e81 648 curv = curvArray[itrack];
649 curverr = curverrArray[itrack];
650 fspTree->Fill();
651 }
652 }
653
367c6d1f 654 if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
f27a7e81 655 curv = 0.5*(curvArray[0]+curvArray[1]);
656 curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
657 fspTree->Fill();
658 }
659 PostData(0,fspTree);
660
a043746c 661 if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return;
f27a7e81 662 // fill ntuple with track-to-track matching
663 Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;
664 Float_t d0[2],z0[2];
665 Float_t sigmad0[2],sigmaz0[2];
666 phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
667 thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
668 phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
669 thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
670 rotymu = TMath::ATan2(track1->Px(),track1->Pz());
671 rotyout = TMath::ATan2(track2->Px(),track2->Pz());
672
673 dphi = phimu - (phiout+TMath::Pi());
674 dtheta = thetamu - (TMath::Pi()-thetaout);
675 if(rotymu>0) {
676 droty = rotymu - (rotyout+TMath::Pi());
677 } else {
678 droty = rotymu - (rotyout-TMath::Pi());
679 }
680
681 Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
682
683 track1->Propagate(alpha,0.,esd->GetMagneticField());
684 track2->Propagate(alpha,0.,esd->GetMagneticField());
685 d0[0] = track1->GetY();
686 z0[0] = track1->GetZ();
687 d0[1] = track2->GetY();
688 z0[1] = track2->GetZ();
689 Float_t dxy = -(d0[0]-d0[1]);
690 Float_t dz = z0[0]-z0[1];
691 sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
692 sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
693 sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
694 sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
695 /*
696 Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
697 track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
698 track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
699 track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
700 track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
701 Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
702 Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
703 Float_t dxaty0 = x1aty0-x2aty0;
704 */
705 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]);
706
707 return;
708}
709
710//________________________________________________________________________
711void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
712{
713 // Extract ITS AliTrackPoints for Cosmics (check angular matching
714 // of top and bottom track, merge the two tracks, if requested)
715 //
716
717 // Set branch addresses for space points tree
718 AliTrackPointArray *arrayForTree=0;
ab6c74ff 719 Float_t curv,curverr,runNumber;
720 TObjString *itsaligndata=0;
3b9267a6 721 TObjString *itscalibrespsdd = 0;
f27a7e81 722 fspTree->SetBranchAddress("SP",&arrayForTree);
723 fspTree->SetBranchAddress("curv",&curv);
724 fspTree->SetBranchAddress("curverr",&curverr);
ab6c74ff 725 fspTree->SetBranchAddress("run",&runNumber);
726 fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
3b9267a6 727 fspTree->SetBranchAddress("ITSCalibRespSDD",&itscalibrespsdd);
ab6c74ff 728
729
730 runNumber = (Float_t)esd->GetRunNumber();
731
732 TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
733 // Get the list of OCDB objects used for reco
734 TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
735 TIter iter2(cdbList);
736 TObjString* cdbEntry=0;
737 TString cdbEntryString;
738 while((cdbEntry =(TObjString*)(iter2.Next()))) {
739 cdbEntryString = cdbEntry->GetString();
740 if(cdbEntryString.Contains("ITS/Align/Data"))
741 itsaligndata = new TObjString(*cdbEntry);
3b9267a6 742 if(cdbEntryString.Contains("ITS/Calib/RespSDD"))
743 itscalibrespsdd = new TObjString(*cdbEntry);
ab6c74ff 744 }
f27a7e81 745
746 Int_t ntracks = esd->GetNumberOfTracks();
747
748 if(ntracks==0) return;
749
750 if(esd->GetPrimaryVertexTracks()->GetNContributors()<=0) return;
751
752 Double_t vtxpos[3]; esd->GetPrimaryVertexTracks()->GetXYZ(vtxpos);
753
754 Int_t ncls=0;
755 Double_t pt=-10000.;
756 Double_t d0z0[2],covd0z0[3];
757 const AliTrackPointArray *array = 0;
758
759 for (Int_t itrack=0; itrack < ntracks; itrack++) {
760 AliESDtrack * track = esd->GetTrack(itrack);
761 if (!track) continue;
762
367c6d1f 763 if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
f27a7e81 764
a043746c 765 if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
766 if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
f27a7e81 767
367c6d1f 768 if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() ||
769 track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
f27a7e81 770
771 pt = track->Pt();
772 ncls = track->GetNcls(0);
773 Double_t maxd=10000.;
774 track->PropagateToDCA(esd->GetPrimaryVertex(),esd->GetMagneticField(),maxd,d0z0,covd0z0);
775
776 // read ITS cluster map
777 Int_t map[6];
778 for(Int_t ilay=0;ilay<6;ilay++) {
779 map[ilay]=0;
780 if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
781 }
782 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()));
783 Int_t idx[12];
784 track->GetITSclusters(idx);
785 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]));
786
787
788 AliTrackPoint point;
789 Int_t ipt,volId,modId,layerId,lay,lad,det;
790 Int_t jpt=0;
791 Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
792
793 array = track->GetTrackPointArray();
794 if(!array) continue;
795 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
796 array->GetPoint(point,ipt);
797 volId = point.GetVolumeID();
798 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
799 if(layerId<1 || layerId>6) continue;
800 if(point.IsExtra() &&
a043746c 801 (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
f27a7e81 802 layerOK[layerId-1]=kTRUE;
803 jpt++;
804 }
805
367c6d1f 806 if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
f27a7e81 807
808 fHistNpoints->Fill(jpt);
809 fHistPt->Fill(pt);
810 PostData(1,fListOfHistos);
811
812 Float_t dzOverlap[2];
813 Double_t globExtra[3],locExtra[3];
814 arrayForTree = new AliTrackPointArray(jpt);
815 jpt=0;
816 array = track->GetTrackPointArray();
817 if(!array) continue;
818 for(ipt=0; ipt<array->GetNPoints(); ipt++) {
819 array->GetPoint(point,ipt);
820 volId = point.GetVolumeID();
821 layerId = AliGeomManager::VolUIDToLayer(volId,modId);
822 if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
823 if(!point.IsExtra() ||
a043746c 824 !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
f27a7e81 825 ncls--;
826 for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
827 AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
828 globExtra[0]=point.GetX();
829 globExtra[1]=point.GetY();
830 globExtra[2]=point.GetZ();
831 AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
832 track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
833 AliTrackPoint pointT;
834 Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
835 for(Int_t lll=0;lll<ipt;lll++) {
836 array->GetPoint(pointT,lll);
837 Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
838 if(layerIdT!=layerId) continue;
839 radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
840 radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
841 phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
842 phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
843 if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
844 thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
845 thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
846 dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
847 if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
848 dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
849 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]);
850 }
851 arrayForTree->AddPoint(jpt,&point);
852 switch(layerId) {
853 case 1:
854 fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
855 break;
856 case 2:
857 fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
858 break;
859 case 3:
860 fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
861 break;
862 case 4:
863 fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
864 break;
865 case 5:
866 fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
867 break;
868 case 6:
869 fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
870 break;
871 }
872 jpt++;
873 }
874
875 curv = track->GetC(esd->GetMagneticField());
876 curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
877
878 fspTree->Fill();
879
880 } // end of tracks loop
881
882 PostData(0,fspTree);
883
884 return;
885}
886
887//________________________________________________________________________
888void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
889{
890 // Terminate analysis
891 //
892 AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
893
894 fspTree = dynamic_cast<TTree*> (GetOutputData(0));
895 if (!fspTree) {
896 printf("ERROR: fspTree not available\n");
897 return;
898 }
899
900 fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
901 if (!fListOfHistos) {
902 printf("ERROR: fListOfHistos not available\n");
903 return;
904 }
905
a043746c 906 fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents"));
f27a7e81 907 fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
908 fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
909 fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
910 fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
911 fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
912 fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
913 fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
914 fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
915 fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
916 fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
917
918
919
920 return;
921}
922