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