ee8703c7166788415c6676d7926cea2f91b45d1e
[u/mrichter/AliRoot.git] / PWG1 / ITS / AliAlignmentDataFilterITS.cxx
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@pd.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 <TMap.h>
31 #include <TVector3.h>
32 #include <TGeoManager.h>
33
34 #include "AliLog.h"
35 #include "AliGeomManager.h"
36 #include "AliITSReconstructor.h"
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
49 ClassImp(AliAlignmentDataFilterITS)
50
51
52 //________________________________________________________________________
53 AliAlignmentDataFilterITS::AliAlignmentDataFilterITS(const char *name):
54 AliAnalysisTask(name,"task"),
55 fOnlySPDFO(kFALSE),
56 fGeometryFileName("geometry.root"),
57 fITSRecoParam(0),
58 fESD(0),
59 fESDfriend(0),
60 fListOfHistos(0),
61 fspTree(0),
62 fHistNevents(0),
63 fHistNpoints(0),
64 fHistPt(0),
65 fHistLayer0(0),
66 fHistLayer1(0),
67 fHistLayer2(0),
68 fHistLayer3(0),
69 fHistLayer4(0),
70 fHistLayer5(0),
71 fntExtra(0),
72 fntCosmicMatching(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 //________________________________________________________________________
86 AliAlignmentDataFilterITS::~AliAlignmentDataFilterITS()
87 {
88   // Destructor
89   if (fListOfHistos) {
90     delete fListOfHistos;
91     fListOfHistos = 0;
92   }
93   if (fspTree) {
94     delete fspTree;
95     fspTree = 0;
96   }
97   if (fHistNevents) {
98     delete fHistNevents;
99     fHistNevents = 0;
100   }
101   if (fHistNpoints) {
102     delete fHistNpoints;
103     fHistNpoints = 0;
104   }
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 }  
142
143 //________________________________________________________________________
144 void 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 {
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     }
190
191     // Disable all branches and enable only the needed ones
192     tree->SetBranchStatus("fTriggerMask", 1);
193     tree->SetBranchStatus("fSPDVertex*", 1);
194     tree->SetBranchStatus("ESDfriend*", 1);
195     tree->SetBranchAddress("ESDfriend.",&fESDfriend);
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 //________________________________________________________________________
211 void AliAlignmentDataFilterITS::Init()
212 {
213   // Initialization
214
215   return;
216 }
217
218 //________________________________________________________________________
219 void AliAlignmentDataFilterITS::CreateOutputObjects()
220 {
221   // Create the output container
222   //
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   }
233
234   // Several histograms are more conveniently managed in a TList
235   fListOfHistos = new TList();
236   fListOfHistos->SetOwner();
237
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);
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
286   fspTree = new TTree("spTree","Tree with ITS track points");
287   const AliTrackPointArray *array = 0;
288   Float_t curv,curverr,runNumber;
289   const TObjString *itsaligndata = 0;
290   fspTree->Branch("SP","AliTrackPointArray",&array);
291   fspTree->Branch("curv",&curv);
292   fspTree->Branch("curverr",&curverr);
293   fspTree->Branch("run",&runNumber);
294   fspTree->Branch("ITSAlignData",&itsaligndata);
295   
296
297   return;
298 }
299
300 //________________________________________________________________________
301 void AliAlignmentDataFilterITS::Exec(Option_t */*option*/)
302 {
303   // Execute analysis for current event:
304   // write ITS AliTrackPoints for selected tracks to fspTree
305   
306   // check the geometry  
307   if(!gGeoManager) { 
308     printf("AliAlignmentDataFilterITS::Exec(): no geometry loaded \n");
309     return;
310   }
311
312   // check if we have AliITSRecoParam
313   if(!GetRecoParam()) {
314     if(!fITSRecoParam) {
315       printf("AliAlignmentDataFilterITS::Exec(): no AliITSRecoParam\n");
316       return;
317     }
318   }
319
320
321   if(!fESD) {
322     printf("AliAlignmentDataFilterITS::Exec(): no ESD \n");
323     return;
324   } 
325   if(!fESDfriend) {
326     printf("AliAlignmentDataFilterITS::Exec(): no ESDfriend \n");
327     return;
328   } 
329   // attach ESDfriend
330   fESD->SetESDfriend(fESDfriend);
331
332   // Post the data for slot 0
333   fHistNevents->Fill(0);
334
335
336   // write field value to spTree UserInfo
337   if(!((fspTree->GetUserInfo())->FindObject("BzkGauss"))) {
338     Double_t bz=fESD->GetMagneticField();
339     TString bzString; bzString+=bz;
340     TObjString *bzObjString = new TObjString(bzString);
341     TList *bzList = new TList();         
342     bzList->SetOwner(1);         
343     bzList->SetName("BzkGauss");         
344     bzList->Add(bzObjString);
345     fspTree->GetUserInfo()->Add(bzList);
346   }
347
348
349   // Process event as Cosmic or Collision
350   //if(esd->GetEventType()== ???? ) {
351   printf("AliAlignmentDataFilterITS::Exec(): MOVE ASAP TO esd->GetEventType() !\n");
352   if(GetRecoParam()->GetAlignFilterCosmics()) {
353     FilterCosmic(fESD);
354   } else {
355     FilterCollision(fESD);
356   }
357
358   PostData(1,fListOfHistos);
359
360   return;
361 }
362
363 //________________________________________________________________________
364 void AliAlignmentDataFilterITS::FilterCosmic(const AliESDEvent *esd)
365 {
366   // Extract ITS AliTrackPoints for Cosmics (check angular matching
367   // of top and bottom track, merge the two tracks, if requested)
368   //
369
370   // Set branch addresses for space points tree
371   AliTrackPointArray *arrayForTree=0;
372   Float_t curv,curverr,runNumber;
373   TObjString *itsaligndata=0;
374   fspTree->SetBranchAddress("SP",&arrayForTree);
375   fspTree->SetBranchAddress("curv",&curv);
376   fspTree->SetBranchAddress("curverr",&curverr);
377   fspTree->SetBranchAddress("run",&runNumber);
378   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
379
380
381   runNumber = (Float_t)esd->GetRunNumber();
382
383   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
384   // Get the list of OCDB objects used for reco 
385   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
386   TIter iter2(cdbList);      
387   TObjString* cdbEntry=0;
388   TString cdbEntryString;
389   while((cdbEntry =(TObjString*)(iter2.Next()))) {
390   cdbEntryString = cdbEntry->GetString();
391   if(cdbEntryString.Contains("ITS/Align/Data")) 
392     itsaligndata = new TObjString(*cdbEntry);
393   }      
394
395
396   TString triggeredClass = fESD->GetFiredTriggerClasses(); 
397   if(fOnlySPDFO && !triggeredClass.Contains("C0SCO-ABCE-NOPF-CENT")) return;
398
399
400   Int_t ntracks = esd->GetNumberOfTracks();
401   if(ntracks<2) return;
402
403   if(esd->GetPrimaryVertexSPD()->GetNContributors()<0) return;
404
405   Double_t vtxpos[3]; esd->GetPrimaryVertexSPD()->GetXYZ(vtxpos);
406
407   Int_t *goodtracksArray = new Int_t[ntracks];
408   Float_t *phiArray = new Float_t[ntracks];
409   Float_t *thetaArray = new Float_t[ntracks];
410   Int_t *nclsArray = new Int_t[ntracks];
411   Int_t ngt=0;
412   Int_t itrack=0;
413   for (itrack=0; itrack < ntracks; itrack++) {
414     AliESDtrack *track = esd->GetTrack(itrack);
415     if (!track) continue;
416
417
418     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
419
420     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
421     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
422
423     Float_t phi = track->GetAlpha()+TMath::ASin(track->GetSnp());
424     Float_t theta = 0.5*TMath::Pi()-TMath::ATan(track->GetTgl());
425
426     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() || 
427        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
428
429     goodtracksArray[ngt] = itrack;
430     phiArray[ngt]        = phi;
431     thetaArray[ngt]      = theta;
432     nclsArray[ngt]       = track->GetNcls(0);
433     ngt++;
434   }
435
436   if(ngt<2) {
437     delete [] goodtracksArray; goodtracksArray=0;
438     delete [] phiArray; phiArray=0;
439     delete [] thetaArray; thetaArray=0;
440     delete [] nclsArray; nclsArray=0;
441     return;
442   }
443
444   // check matching of the two tracks from the muon
445   Float_t min = 10000000.;
446   Int_t maxCls = 0;
447   Int_t good1 = -1, good2 = -1;
448   for(Int_t itr1=0; itr1<ngt-1; itr1++) {
449     for(Int_t itr2=itr1+1; itr2<ngt; itr2++) {
450       Float_t deltatheta = TMath::Abs(TMath::Pi()-thetaArray[itr1]-thetaArray[itr2]);
451       if(deltatheta>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
452       Float_t deltaphi = TMath::Abs(TMath::Abs(phiArray[itr1]-phiArray[itr2])-TMath::Pi());
453       if(deltaphi>GetRecoParam()->GetAlignFilterMaxMatchingAngle()) continue;
454       if(nclsArray[itr1]+nclsArray[itr2] > maxCls) {
455         maxCls = nclsArray[itr1]+nclsArray[itr2];
456         min = deltatheta+deltaphi;
457         good1 = goodtracksArray[itr1];
458         good2 = goodtracksArray[itr2];
459       } else if(nclsArray[itr1]+nclsArray[itr2] == maxCls) {
460         if(deltatheta+deltaphi < min) {
461           min = deltatheta+deltaphi;
462           good1 = goodtracksArray[itr1];
463           good2 = goodtracksArray[itr2];
464         }
465       }
466     }
467   }
468   
469   delete [] goodtracksArray; goodtracksArray=0;
470   delete [] phiArray; phiArray=0;
471   delete [] thetaArray; thetaArray=0;
472   delete [] nclsArray; nclsArray=0;
473
474   if(good1<0) return;
475   AliDebug(2,"ok track matching");
476   
477   // track1 will be the inward track (top)
478   // track2 the outward (bottom)
479   AliESDtrack *track1=0; 
480   AliESDtrack *track2=0;
481   AliESDtrack *track = esd->GetTrack(good1);
482   if(track->Py()>0) { 
483     track1 = esd->GetTrack(good1);
484     track2 = esd->GetTrack(good2);
485   } else {
486     track1 = esd->GetTrack(good2);
487     track2 = esd->GetTrack(good1);
488   }
489
490   AliTrackPoint point;
491   const AliTrackPointArray *array=0;
492   Int_t ipt,volId,modId,layerId,lay,lad,det;
493   Int_t jpt=0;
494   Bool_t layerOK[6][2]; 
495   Int_t nclsTrk[2]={0,0};
496
497   for(Int_t l1=0;l1<6;l1++) for(Int_t l2=0;l2<2;l2++) layerOK[l1][l2]=kFALSE;
498     
499   for(itrack=0; itrack<2; itrack++) {
500     if(itrack==0) {
501       track = track1;
502     } else {
503       track = track2;
504     }
505     array = track->GetTrackPointArray();
506     if(!array) {
507       AliWarning("No tracks points avaialble");
508       continue;
509     }
510     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
511       array->GetPoint(point,ipt);
512       volId = point.GetVolumeID();
513       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
514       AliDebug(2,Form("%d %d\n",ipt,layerId-1));
515       if(point.IsExtra() && 
516          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
517       if(layerId>6) continue;
518       if(!GetRecoParam()->GetAlignFilterUseLayer(layerId-1)) continue;
519       // check minAngleWrtITSModulePlanes
520       Double_t p[3]; track->GetDirection(p);
521       TVector3 pvec(p[0],p[1],p[2]);
522       Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
523       TVector3 normvec(rot[1],rot[4],rot[7]);
524       Double_t angle = pvec.Angle(normvec);
525       if(angle>0.5*TMath::Pi()) angle = TMath::Pi()-angle;
526       angle = 0.5*TMath::Pi()-angle;
527       if(angle<GetRecoParam()->GetAlignFilterMinAngleWrtModulePlanes()) continue;
528       layerOK[layerId-1][itrack]=kTRUE;
529       jpt++;
530       nclsTrk[itrack]++;
531     }
532   }
533   AliDebug(2,Form("nClsTrk1 %d nClsTrk2 %d\n",nclsTrk[0],nclsTrk[1]));
534     
535   // read ITS cluster maps
536   Int_t map1[6],map2[6];
537   for(Int_t ilay=0;ilay<6;ilay++) {
538     map1[ilay]=0; map2[ilay]=0;
539     if(track1->HasPointOnITSLayer(ilay)) map1[ilay]=1;
540     if(track2->HasPointOnITSLayer(ilay)) map2[ilay]=1;
541   }
542   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()));
543   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()));
544   Int_t idx1[12],idx2[12];
545   track1->GetITSclusters(idx1);
546   track2->GetITSclusters(idx2);
547   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]));
548   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]));
549   
550
551   if(jpt<GetRecoParam()->GetAlignFilterMinITSPointsMerged()) return;
552   AliDebug(2,Form(" Total points %d, accepted\n",jpt));  
553   fHistNpoints->Fill(jpt);
554   fHistPt->Fill(0.5*(track1->Pt()+track2->Pt()));
555   
556   Float_t d0z0mu[2];
557   track1->GetDZ(0,0,0,esd->GetMagneticField(),d0z0mu);
558   //printf("d0mu %f  z0mu %f\n",d0z0mu[0],d0z0mu[1]);
559
560   Float_t dzOverlap[2];
561   Float_t curvArray[2],curverrArray[2];
562   Double_t globExtra[3],locExtra[3];
563   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) 
564     arrayForTree = new AliTrackPointArray(jpt);
565   
566   jpt=0;
567   for(itrack=0; itrack<2; itrack++) {
568     if(itrack==0) {
569       track = track1;
570     } else {
571       track = track2;
572     }
573     curvArray[itrack] = track->GetC(esd->GetMagneticField());
574     curverrArray[itrack] = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
575
576     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
577       jpt=0;
578       arrayForTree = new AliTrackPointArray(nclsTrk[itrack]);
579     }
580     array = track->GetTrackPointArray();
581     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
582       array->GetPoint(point,ipt);
583       volId = point.GetVolumeID();
584       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
585       if(layerId>6 || !layerOK[layerId-1][itrack]) continue;
586       arrayForTree->AddPoint(jpt,&point);
587       jpt++;
588       switch(layerId) {
589       case 1:
590         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
591         break;
592       case 2:
593         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
594         break;
595       case 3:
596         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
597         break;
598       case 4:
599         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
600         break;
601       case 5:
602         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
603         break;
604       case 6:
605         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
606         break;
607       }
608       // Post the data for slot 0
609       if(jpt==1) PostData(1,fListOfHistos); // only if this is the first points
610       if(!point.IsExtra() || 
611          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
612       nclsTrk[itrack]--;
613       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
614       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
615       globExtra[0]=point.GetX();
616       globExtra[1]=point.GetY();
617       globExtra[2]=point.GetZ();
618       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
619       //printf("%d %d %d %d %d  %f %f %f\n",volId,modId,lay,lad,det,locExtra[0],locExtra[1],locExtra[2]);
620       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
621       AliTrackPoint pointT;
622       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
623       for(Int_t lll=0;lll<ipt;lll++) {
624         array->GetPoint(pointT,lll);
625         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
626         if(layerIdT!=layerId) continue;
627         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
628         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
629         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
630         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
631         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
632         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
633         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
634         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
635         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
636         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
637         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]);
638       }
639     }
640
641     if(!(GetRecoParam()->GetAlignFilterCosmicMergeTracks())) {
642       curv = curvArray[itrack];
643       curverr = curverrArray[itrack];
644       fspTree->Fill();
645     }
646   }
647
648   if(GetRecoParam()->GetAlignFilterCosmicMergeTracks()) {
649     curv = 0.5*(curvArray[0]+curvArray[1]);
650     curverr = 0.5*TMath::Sqrt(curverrArray[0]*curverrArray[0]+curverrArray[1]*curverrArray[1]);
651     fspTree->Fill();
652   }
653   PostData(0,fspTree);
654
655   if(!(GetRecoParam()->GetAlignFilterFillQANtuples())) return; 
656   // fill ntuple with track-to-track matching
657   Float_t phimu,thetamu,phiout,thetaout,dphi,dtheta,rotymu,rotyout,droty;    
658   Float_t d0[2],z0[2];
659   Float_t sigmad0[2],sigmaz0[2];
660   phimu = track1->GetAlpha()+TMath::ASin(track1->GetSnp());
661   thetamu = 0.5*TMath::Pi()-TMath::ATan(track1->GetTgl());
662   phiout = track2->GetAlpha()+TMath::ASin(track2->GetSnp());
663   thetaout = 0.5*TMath::Pi()-TMath::ATan(track2->GetTgl());
664   rotymu = TMath::ATan2(track1->Px(),track1->Pz());
665   rotyout = TMath::ATan2(track2->Px(),track2->Pz());
666
667   dphi = phimu - (phiout+TMath::Pi());
668   dtheta = thetamu - (TMath::Pi()-thetaout);
669   if(rotymu>0) {
670     droty = rotymu - (rotyout+TMath::Pi());
671   } else {
672     droty = rotymu - (rotyout-TMath::Pi());
673   }
674
675   Double_t alpha = TMath::ATan2(track1->Py(),track1->Px());
676
677   track1->Propagate(alpha,0.,esd->GetMagneticField());
678   track2->Propagate(alpha,0.,esd->GetMagneticField());
679   d0[0] = track1->GetY();
680   z0[0] = track1->GetZ();
681   d0[1] = track2->GetY();
682   z0[1] = track2->GetZ();
683   Float_t dxy = -(d0[0]-d0[1]);
684   Float_t dz  = z0[0]-z0[1];
685   sigmad0[0] = TMath::Sqrt(track1->GetSigmaY2());
686   sigmaz0[0] = TMath::Sqrt(track1->GetSigmaZ2());
687   sigmad0[1] = TMath::Sqrt(track2->GetSigmaY2());
688   sigmaz0[1] = TMath::Sqrt(track2->GetSigmaZ2());
689   /*  
690   Double_t xyz1atxl0[3],xyz1atxl1[3],xyz2atxl0[3],xyz2atxl1[3];
691   track1->GetXYZAt(0.,esd->GetMagneticField(),xyz1atxl0);
692   track1->GetXYZAt(1.,esd->GetMagneticField(),xyz1atxl1);
693   track2->GetXYZAt(0.,esd->GetMagneticField(),xyz2atxl0);
694   track2->GetXYZAt(1.,esd->GetMagneticField(),xyz2atxl1);
695   Float_t x1aty0 = (xyz1atxl0[0]*xyz1atxl1[1]-xyz1atxl0[1]*xyz1atxl1[0])/(xyz1atxl1[1]-xyz1atxl0[1]);
696   Float_t x2aty0 = (xyz2atxl0[0]*xyz2atxl1[1]-xyz2atxl0[1]*xyz2atxl1[0])/(xyz2atxl1[1]-xyz2atxl0[1]);
697   Float_t dxaty0 = x1aty0-x2aty0;
698   */
699   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]);
700   
701   return;
702 }
703
704 //________________________________________________________________________
705 void AliAlignmentDataFilterITS::FilterCollision(const AliESDEvent *esd)
706 {
707   // Extract ITS AliTrackPoints for Cosmics (check angular matching
708   // of top and bottom track, merge the two tracks, if requested)
709   //
710
711   // Set branch addresses for space points tree
712   AliTrackPointArray *arrayForTree=0;
713   Float_t curv,curverr,runNumber;
714   TObjString *itsaligndata=0;
715   fspTree->SetBranchAddress("SP",&arrayForTree);
716   fspTree->SetBranchAddress("curv",&curv);
717   fspTree->SetBranchAddress("curverr",&curverr);
718   fspTree->SetBranchAddress("run",&runNumber);
719   fspTree->SetBranchAddress("ITSAlignData",&itsaligndata);
720
721
722   runNumber = (Float_t)esd->GetRunNumber();
723
724   TTree* esdTree = dynamic_cast<TTree*> (GetInputData(0));
725   // Get the list of OCDB objects used for reco 
726   TList *cdbList = (TList*)(esdTree->GetTree()->GetUserInfo())->FindObject("cdbList");
727   TIter iter2(cdbList);      
728   TObjString* cdbEntry=0;
729   TString cdbEntryString;
730   while((cdbEntry =(TObjString*)(iter2.Next()))) {
731   cdbEntryString = cdbEntry->GetString();
732   if(cdbEntryString.Contains("ITS/Align/Data")) 
733     itsaligndata = new TObjString(*cdbEntry);
734   }      
735
736   Int_t ntracks = esd->GetNumberOfTracks();
737
738   if(ntracks==0) return;
739
740   if(esd->GetPrimaryVertexTracks()->GetNContributors()<=0) return;
741
742   Double_t vtxpos[3]; esd->GetPrimaryVertexTracks()->GetXYZ(vtxpos);
743
744   Int_t ncls=0;
745   Double_t pt=-10000.;
746   Double_t d0z0[2],covd0z0[3];
747   const AliTrackPointArray *array = 0;
748
749   for (Int_t itrack=0; itrack < ntracks; itrack++) {
750     AliESDtrack * track = esd->GetTrack(itrack);
751     if (!track) continue;
752
753     if(track->GetNcls(0)<GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
754
755     if((GetRecoParam()->GetAlignFilterOnlyITSSATracks()) && track->GetNcls(1)>0) continue;
756     if((GetRecoParam()->GetAlignFilterOnlyITSTPCTracks()) && track->GetNcls(1)==0) continue;
757
758     if(track->Pt()<GetRecoParam()->GetAlignFilterMinPt() || 
759        track->Pt()>GetRecoParam()->GetAlignFilterMaxPt()) continue;
760
761     pt = track->Pt();
762     ncls = track->GetNcls(0);
763     Double_t maxd=10000.;
764     track->PropagateToDCA(esd->GetPrimaryVertex(),esd->GetMagneticField(),maxd,d0z0,covd0z0);
765
766     // read ITS cluster map
767     Int_t map[6];
768     for(Int_t ilay=0;ilay<6;ilay++) {
769       map[ilay]=0;
770       if(track->HasPointOnITSLayer(ilay)) map[ilay]=1;
771     }
772     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()));
773     Int_t idx[12];
774     track->GetITSclusters(idx);
775     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]));
776     
777
778     AliTrackPoint point;
779     Int_t ipt,volId,modId,layerId,lay,lad,det;
780     Int_t jpt=0;
781     Bool_t layerOK[6]; for(Int_t l1=0;l1<6;l1++) layerOK[l1]=kFALSE;
782     
783     array = track->GetTrackPointArray();
784     if(!array) continue;
785     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
786       array->GetPoint(point,ipt);
787       volId = point.GetVolumeID();
788       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
789       if(layerId<1 || layerId>6) continue;
790       if(point.IsExtra() && 
791          (GetRecoParam()->GetAlignFilterSkipExtra())) continue;
792       layerOK[layerId-1]=kTRUE;
793       jpt++;
794     }
795
796     if(jpt < GetRecoParam()->GetAlignFilterMinITSPoints()) continue;
797
798     fHistNpoints->Fill(jpt);
799     fHistPt->Fill(pt);
800     PostData(1,fListOfHistos);
801
802     Float_t dzOverlap[2];
803     Double_t globExtra[3],locExtra[3];
804     arrayForTree = new AliTrackPointArray(jpt);
805     jpt=0;
806     array = track->GetTrackPointArray();
807     if(!array) continue;
808     for(ipt=0; ipt<array->GetNPoints(); ipt++) {
809       array->GetPoint(point,ipt);
810       volId = point.GetVolumeID();
811       layerId = AliGeomManager::VolUIDToLayer(volId,modId);
812       if(layerId<1 || layerId>6 || !layerOK[layerId-1]) continue;
813       if(!point.IsExtra() || 
814          !(GetRecoParam()->GetAlignFilterFillQANtuples())) continue;
815       ncls--;
816       for(Int_t ll=1;ll<layerId;ll++) modId+=AliITSgeomTGeo::GetNLadders(ll)*AliITSgeomTGeo::GetNDetectors(ll);
817       AliITSgeomTGeo::GetModuleId(modId,lay,lad,det);
818       globExtra[0]=point.GetX();
819       globExtra[1]=point.GetY();
820       globExtra[2]=point.GetZ();
821       AliITSgeomTGeo::GlobalToLocal(lay,lad,det,globExtra,locExtra);
822       track->GetDZ(point.GetX(),point.GetY(),point.GetZ(),esd->GetMagneticField(),dzOverlap);
823       AliTrackPoint pointT;
824       Float_t radius,radiusT,phiv,phivT,thetav,thetavT;
825       for(Int_t lll=0;lll<ipt;lll++) {
826         array->GetPoint(pointT,lll);
827         Int_t layerIdT = AliGeomManager::VolUIDToLayer(pointT.GetVolumeID(),modId);
828         if(layerIdT!=layerId) continue;
829         radius=TMath::Sqrt((point.GetX()-vtxpos[0])*(point.GetX()-vtxpos[0])+(point.GetY()-vtxpos[1])*(point.GetY()-vtxpos[1]));
830         radiusT=TMath::Sqrt((pointT.GetX()-vtxpos[0])*(pointT.GetX()-vtxpos[0])+(pointT.GetY()-vtxpos[1])*(pointT.GetY()-vtxpos[1]));
831         phiv=TMath::ATan2(point.GetY()-vtxpos[1],point.GetX()-vtxpos[0]);
832         phivT=TMath::ATan2(pointT.GetY()-vtxpos[1],pointT.GetX()-vtxpos[0]);
833         if(TMath::Abs(point.GetZ()-vtxpos[2])<0.00001 || TMath::Abs(pointT.GetZ()-vtxpos[2])<0.00001) continue;
834         thetav=TMath::ATan(radius/(point.GetZ()-vtxpos[2]));
835         thetavT=TMath::ATan(radiusT/(pointT.GetZ()-vtxpos[2]));
836         dzOverlap[0]=(Float_t)((phivT-phiv)*0.5*(radiusT+radius));
837         if(TMath::Abs(TMath::Tan(0.5*(thetav+thetavT)))<0.00001) continue;
838         dzOverlap[1]=(Float_t)((pointT.GetZ()-point.GetZ())-(radiusT-radius)/TMath::Tan(0.5*(thetav+thetavT)));
839         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]);
840       }
841       arrayForTree->AddPoint(jpt,&point);
842       switch(layerId) {
843       case 1:
844         fHistLayer0->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
845         break;
846       case 2:
847         fHistLayer1->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
848         break;
849       case 3:
850         fHistLayer2->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
851         break;
852       case 4:
853         fHistLayer3->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
854         break;
855       case 5:
856         fHistLayer4->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
857         break;
858       case 6:
859         fHistLayer5->Fill(TMath::ATan2(point.GetY(),point.GetX()),point.GetZ());
860         break;
861       }
862       jpt++;
863     }
864
865     curv = track->GetC(esd->GetMagneticField());
866     curverr = TMath::Sqrt(track->GetSigma1Pt2())*track->GetC(esd->GetMagneticField())/track->OneOverPt();
867
868     fspTree->Fill();
869  
870   } // end of tracks loop
871
872   PostData(0,fspTree);
873
874   return;
875 }
876
877 //________________________________________________________________________
878 void AliAlignmentDataFilterITS::Terminate(Option_t */*option*/)
879 {
880   // Terminate analysis
881   //
882   AliDebug(2,"AliITSAlignmentDataFiler: Terminate() \n");
883
884   fspTree = dynamic_cast<TTree*> (GetOutputData(0));
885   if (!fspTree) {     
886     printf("ERROR: fspTree not available\n");
887     return;
888   }
889
890   fListOfHistos = dynamic_cast<TList*> (GetOutputData(1));
891   if (!fListOfHistos) {     
892     printf("ERROR: fListOfHistos not available\n");
893     return;
894   }
895
896   fHistNevents = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNevents"));
897   fHistNpoints = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistNpoints"));
898   fHistPt = dynamic_cast<TH1F*>(fListOfHistos->FindObject("fHistPt"));
899   fHistLayer0 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer0"));
900   fHistLayer1 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer1"));
901   fHistLayer2 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer2"));
902   fHistLayer3 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer3"));
903   fHistLayer4 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer4"));
904   fHistLayer5 = dynamic_cast<TH2F*>(fListOfHistos->FindObject("fHistLayer5"));
905   fntExtra = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntExtra"));
906   fntCosmicMatching = dynamic_cast<TNtuple*>(fListOfHistos->FindObject("fntCosmicMatching"));
907
908
909
910   return;
911 }
912