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