eac8e3abbdf65cfa5700cce8eae4b522002d4fe8
[u/mrichter/AliRoot.git] / ITS / UPGRADE / v0 / residuals.C
1
2
3
4 Int_t nLayers =0;
5 Int_t GetLayer(AliTrackReference *ref, TObject *obj);
6
7 void residuals(){
8
9   gSystem->Load("libITSUpgradeBase.so");
10   gSystem->Load("libITSUpgradeRec.so");
11   gSystem->Load("libITSUpgradeSim.so");
12
13   gROOT->SetStyle("Plain");
14   gStyle->SetOptStat(1111111);
15   Int_t nbins=100;
16   Int_t xmin=0;
17   Int_t xmax=50000;//00*1e-09;
18
19   AliITSsegmentationUpgrade *seg = new AliITSsegmentationUpgrade();
20   
21   if(!seg){
22     printf("no segmentation info available... Exiting");
23     return;
24   }
25
26   nLayers = seg->GetNLayers();
27
28   Int_t nbins=100;
29   Int_t xmin = -10;
30   Int_t xmax= 10;
31   Int_t zmin = -3000;
32   Int_t zmax= 3000;
33
34   TH1D **hDiffx, **hDiffy, **hDiffz, **hDiffRphi, **hDiffR;
35   
36   hDiffx    = new TH1D[nLayers];
37   hDiffy    = new TH1D[nLayers];
38   hDiffz    = new TH1D[nLayers];
39   hDiffRphi = new TH1D[nLayers];
40   hDiffR    = new TH1D[nLayers];
41     
42   for(Int_t i=0; i<nLayers; i++){
43     hDiffx[i] = new TH1D(Form("hDiffX%i",i),Form("#Delta X Glob (Clusters - track ref  [Layer %i])",i),100,-40,40);
44     hDiffy[i] = new TH1D(Form("hDiffY%i",i),Form("#Delta Y Glob (Clusters - track ref  [Layer %i])",i),100,-40,40);
45     if(i<4) hDiffz[i] = new TH1D(Form("hDiffZ%i",i),Form("#Delta Z Glob (Clusters - track ref  [Layer %i])",i),50,-300,300);
46     else hDiffz[i] = new TH1D(Form("hDiffZ%i",i),Form("#Delta Z Glob (Clusters - track ref  [Layer %i])",i),100,-2000,2000);
47     hDiffRphi[i] = new TH1D(Form("hDiffRphi%i",i),Form("#Delta R #phi [Layer %i]",i),200,-100,100);
48     hDiffRphi[i]->SetXTitle("#mum");
49     hDiffR[i]= new TH1D(Form("hR%i",i),Form("#Delta R (cluster-track ref) [Layer %i]",i),nbins,xmin/20.,xmax/20.);
50     hDiffR[i]->SetXTitle("#mum");
51   }
52
53   gAlice=NULL;
54   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
55   runLoader->LoadgAlice();
56
57   gAlice = runLoader->GetAliRun();
58
59   runLoader->LoadHeader();
60   runLoader->LoadKinematics();
61   runLoader->LoadTrackRefs();
62   runLoader->LoadRecPoints();
63
64   AliLoader *dl = runLoader->GetDetectorLoader("ITS");
65
66   //Trackf
67   TTree *trackRefTree = 0x0; 
68   TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000);
69   //RECPOINT
70   TTree *clusTree = 0x0;
71   TClonesArray statITSCluster("AliITSRecPointU");
72
73
74   AliITSsegmentationUpgrade *segmentation2=segmentation2 = new AliITSsegmentationUpgrade();
75   // Loop over events
76   for(Int_t i=0; i<runLoader->GetNumberOfEvents(); i++){
77     runLoader->GetEvent(i); 
78     AliStack *stack = runLoader->Stack(); 
79     trackRefTree=runLoader->TreeTR();
80     TBranch *br = trackRefTree->GetBranch("TrackReferences");
81     if(!br) {
82       printf("no TR branch available , exiting \n");
83       return;
84     }
85     br->SetAddress(&trackRef);
86
87     printf("event : %i \n",i);
88     // init the clusters tree 
89     clusTree=dl->TreeR();
90     TClonesArray *ITSCluster = &statITSCluster;
91     TBranch* itsClusterBranch=clusTree->GetBranch("ITSRecPoints");
92     if (!itsClusterBranch) {
93       printf("can't get the branch with the ITS clusters ! \n");
94       return;
95     }
96
97     itsClusterBranch->SetAddress(&ITSCluster);
98
99     // init the trackRef tree 
100     trackRefTree=runLoader->TreeTR();
101     trackRefTree->SetBranchAddress("TrackReferences",&trackRef);
102
103     //for(Int_t il=0; il<nLayers; il++){
104     if (!clusTree->GetEvent(0))    continue;
105     Int_t nCluster = ITSCluster->GetEntriesFast();
106     for(Int_t in=0; in<nCluster; in++){
107       AliITSRecPointU *recp = (AliITSRecPointU*)ITSCluster->UncheckedAt(in);
108       Int_t il = recp->GetLayer();  
109       Double_t xr,yr,zr = 0.;
110       Double_t xz[2];
111       xz[0]= recp->GetDetLocalX(); //gets fXloc
112       xz[1]= recp->GetDetLocalZ();   //gets fZloc
113       segmentation2->DetToGlobal(il,recp->GetModule(),xz[0], xz[1],xr,yr,zr);
114       for(Int_t iLabel =0; iLabel<12; iLabel++){
115         if(recp->GetTrackID(iLabel)<0) continue; 
116         trackRefTree->GetEntry(stack->TreeKEntry(recp->GetTrackID(iLabel)));    
117         Int_t nref=trackRef->GetEntriesFast();
118         for(Int_t iref =0; iref<nref; iref++){
119           Double_t x,y,z=0.;
120           Int_t labelTR=-999;
121           AliTrackReference *trR = (AliTrackReference*)trackRef->At(iref);
122           if(!trR) continue;
123           if(trR->DetectorId()!=AliTrackReference::kITS) continue;
124           if(!stack->IsPhysicalPrimary(recp->GetTrackID(iLabel))) continue;   
125           Int_t lay = GetLayer(trR, seg);    
126           if(TMath::Abs(lay-il)>0.5) continue;          
127           x=trR->X();
128           y=trR->Y();
129           z=trR->Z();
130           Double_t xx, zz;
131           Int_t mod;
132           segmentation2->GlobalToDet(il,x,y,x,xx,zz,mod);
133           //printf("cluster module %i  - module %i \n",recp->GetModule(),mod);
134           Double_t rTr = TMath::Sqrt(x*x+y*y); 
135
136           Double_t difx, dify, difz=0.;
137           Double_t dify=0.;
138           Double_t difz=0.;
139           Double_t deltaR=0.;
140           difx = xr - x;
141           dify = yr - y;
142           difz = zr - z;
143           Double_t phiGlob = TMath::ATan2(yr,xr);
144           if(yr<0) phiGlob+=TMath::TwoPi();
145           Double_t phiTr = TMath::ATan2(y,x);
146           if(y<0) phiTr+=TMath::TwoPi();
147           Double_t deltaPhi =  phiGlob-phiTr;
148               
149           Double_t deltaRphi = (phiGlob*TMath::Sqrt(xr*xr+yr*yr)-phiTr*rTr)*1e+04;
150           if(TMath::Abs(deltaRphi)/1e+04 > 3*(10000.*seg->GetCellSizeX(lay))) {
151             // printf("   Layer %i   type %i  dRphi %f  segm %f  \n",lay,recp->GetType(),deltaRphi/1e+04,seg->GetCellSizeX(lay)*10000.);
152             // printf("   Layer %i   type %i   Track (x,y,z) = (%f,%f,%f) Phi %f   Rphi %f         Cluster (x,y,z) = (%f,%f,%f)  Phi %f Rphi %f \n",lay,recp->GetType(),x,y,z,TMath::RadToDeg()*phiTr,phiTr*TMath::Sqrt(x*x+y*y),xr,yr,zr,TMath::RadToDeg()*phiGlob,phiGlob*TMath::Sqrt(xr*xr+yr*yr));
153           }              
154           deltaR = (TMath::Sqrt(xr*xr+yr*yr)-TMath::Sqrt(x*x+y*y))*1e+04;
155
156           hDiffx[il]->Fill(difx*1e04);
157           hDiffy[il]->Fill(dify*1e04);
158           hDiffz[il]->Fill(difz*1e04);
159           hDiffRphi[il]->Fill(deltaRphi);
160           hDiffR[il]->Fill(deltaR);
161            
162         }//loop entries fast 
163       } 
164     }//loop clusters
165     //}//loop layer
166   }//loop entries runloader
167   
168   TCanvas *c3x = new TCanvas("cGx","Delta in X global",200,10,900,900);
169   c3x->Divide(3,nLayers/3);
170   TCanvas *c3y = new TCanvas("cGy","Delta in Y global",200,10,900,900);
171   c3y->Divide(3,nLayers/3);
172   TCanvas *c3z = new TCanvas("cGz","Delta in Z global",200,10,900,900);
173   c3z->Divide(3,nLayers/3);
174   TCanvas *c3rphi = new TCanvas("cRphi","Delta R-Phi",200,10,900,900);
175   c3rphi->Divide(3,nLayers/3);
176   TCanvas *c3r = new TCanvas("cR","Delta R-Phi",200,10,900,900);
177   c3r->Divide(3,nLayers/3);
178   
179
180   for(Int_t iL=0; iL< nLayers; iL++){
181     c3x->cd(iL+1);
182     hDiffx[iL]->Draw();
183     c3y->cd(iL+1);
184     hDiffy[iL]->Draw();
185     c3z->cd(iL+1);
186     hDiffz[iL]->Draw();
187     c3rphi->cd(iL+1);
188     hDiffRphi[iL]->Draw();
189     c3r->cd(iL+1);
190     hDiffR[iL]->Draw();
191   }
192  
193 }// main macro read
194
195
196 Int_t GetLayer(AliTrackReference *ref, TObject *obj){
197  
198   AliITSsegmentationUpgrade *seg = (AliITSsegmentationUpgrade*)obj;
199   Int_t ilayer =-1;
200   for(Int_t ila=0; ila<6; ila++){
201     if(ila==(nLayers-1)){
202       if(ref->R()>seg->GetRadius(ila)+seg->GetThickness(ila)) ilayer=(nLayers-1);
203     }
204     else {
205       if(ref->R()>seg->GetRadius(ila)+seg->GetThickness(ila) && ref->R()<seg->GetRadius(ila+1)) ilayer=ila;
206     }
207   }
208   return ilayer;
209 }
210