store also difference in local Y
[u/mrichter/AliRoot.git] / ITS / RealignmentAnalysisITSIndependVol.C
1 #include <TArray.h>
2 #include <TFile.h>
3 #include <TStopwatch.h>
4 #include <TArray.h>
5 #include <TNtuple.h>
6 #include <TGeoManager.h>
7 #include "AliGeomManager.h"
8 #include "AliAlignmentTracks.h"
9 #include "AliTrackFitter.h"
10 #include "AliTrackFitterKalman.h"
11 #include "AliTrackFitterRieman.h"
12 #include "AliTrackResidualsFast.h"
13 #include "AliTrackResidualsChi2.h"
14 #include "AliTrackResidualsLinear.h"
15
16
17 void RealignmentITSIndependVol(TString minimizer="fast",const Int_t fit=0,const Int_t iter1=1,const Int_t iterations=5,const Int_t minNtracks=-10,const Int_t layer=0,const Int_t minTrackPoint=6,TString fileintro="AliTrackPoints.root",TString geometryfile="geometry.root"){
18   
19   // minimizer="fast"->AliTrackResidualFast minimizer
20   //           "minuit"->AliTrackResidualChi2 minimizer
21   //           "minuitnorot"->AliTrackResidualChi2 minimizer without rotations degrees of freedom
22   //           "linear"->AliTrackResidualLinear minimizer    
23   //fit=0-> Riemann Fitter, fit=1->Kalman
24   //iter1=#iterations inside AliAlignmentTracks::AliAlignVolumes method 
25   //iterations=#iterations of the entire procedure
26   //layer=0->all ITS, otherways the usual notation is considered (1=SPD1,2=SPD2,3=SDD1,4=SDD2,5=SSD1,6=SSD2)
27   //minNtracks=minimun number of tracks passing through a module in order to try to realign the module itsself
28   //           if minNtracks<0, minimun number of tracks is |minNtracks|*minNumPoint[layer]/fact (see the code below): this allows a different
29   //           choice of the number of tracks required on different layers and to vary these numbers once tuned the relative proportions.  
30   //minTrackPoint=minimun number of "good" points required to a track (THE POINT ON THE MODULE THAT IS GOING TO BE REALIGNED 
31   //              IS NEVER CONSIDERED->max number can be required is 11 for cosmics tracks) for the track being considered in the minimization
32   //fileintro=file into which the Tree with the space points is stored
33   //geometryfile=file containing the geometry
34
35
36   TArrayI volIDs2(2200); 
37   volIDs2.Reset(0);
38   TArrayI volIDs(1);
39   TString outname;
40   Int_t layerNum,modNum,iLayer,iLayerToAlign,j=0,size=0,lastVolid=0;
41   Int_t minNumPoint[6]={200,200,200,200,100,100}; 
42   Double_t fact=10; 
43   TNtuple *ntVolumeAlign=new TNtuple("ntVolumeAlign","NTuple with volume tried to be realigned","layerNum:modNum:volumeIDnum"); 
44
45
46   AliAlignmentTracks *AliAlTrack=new AliAlignmentTracks();
47  
48   AliAlTrack->SetPointsFilename(fileintro.Data());
49
50   AliTrackFitter *fitter;
51   if(fit==1)fitter= new AliTrackFitterKalman();
52   else fitter=new AliTrackFitterRieman();
53
54   fitter->SetMinNPoints(minTrackPoint);
55
56   AliAlTrack->SetTrackFitter(fitter);
57
58
59   AliTrackResiduals *res;
60   if(minimizer=="minuit"){
61     res = new AliTrackResidualsChi2();
62   }
63   else if(minimizer=="minuitnorot"){
64     res = new AliTrackResidualsChi2();
65     res->FixParameter(3);
66     res->FixParameter(4);
67     res->FixParameter(5);
68   }
69
70   else if(minimizer=="fast"){
71     res = new AliTrackResidualsFast();
72   } 
73   else if(minimizer=="linear"){
74     res = new AliTrackResidualsLinear();
75   } 
76
77   else {
78     printf("Trying to set a non existing minimizer! \n");
79     return;
80   }
81
82   res->SetMinNPoints(1);
83   AliAlTrack->SetMinimizer(res);
84   
85   if(!gGeoManager) TGeoManager::Import(geometryfile.Data());
86    
87   
88   TStopwatch *timer=new TStopwatch();
89   timer->Start(); 
90   AliAlTrack->BuildIndex();
91   j=0;
92   UShort_t volid;
93   for(Int_t iter=0;iter<iterations;iter++){
94     for(iLayerToAlign=1;iLayerToAlign<=6;iLayerToAlign++){
95       if(layer!=0&&(Int_t)iLayerToAlign!=layer)continue;
96       j=0;
97       size=0;
98       for(Int_t k=1;k<=6;k++){
99         size+=(Int_t)AliGeomManager::LayerSize(k);
100         printf("size: %d \n",size);
101       }
102   
103       
104       for (Int_t iModule=0;iModule<AliGeomManager::LayerSize(iLayerToAlign);iModule++){      
105         j=0;
106         if(minNtracks<0){       
107           if(AliAlTrack->GetLastIndex(iLayerToAlign-1,iModule)<minNumPoint[iLayerToAlign-1]*(-1*minNtracks/fact))continue;
108         }       
109         else if(AliAlTrack->GetLastIndex(iLayerToAlign-1,iModule)<minNtracks)continue;
110
111         UShort_t volidAl = AliGeomManager::LayerToVolUID(iLayerToAlign,iModule);
112         
113         TArrayI volIDsFit(size-1);
114         for (iLayer=(Int_t)AliGeomManager::kSPD1;iLayer<(Int_t)AliGeomManager::kTPC1;iLayer++){
115           for (Int_t iModule2=0;iModule2<AliGeomManager::LayerSize(iLayer);iModule2++){
116             volid = AliGeomManager::LayerToVolUID(iLayer,iModule2);
117             if(volid==volidAl)continue;
118             volIDsFit.AddAt(volid,j);
119             j++;
120           }
121         }
122         volIDs.AddAt((Int_t)volidAl,0);
123         if(iter==iterations-1){
124           volIDs2.AddAt(volidAl,lastVolid);
125           lastVolid++;
126         }
127         AliAlTrack->AlignVolumes(&volIDs,&volIDsFit,AliGeomManager::kSPD1,AliGeomManager::kSSD2,iter1);  
128       }
129     }
130     
131     if((iter+1)%5==0||iter==0||iter==1||iter==2||iter==3||iter==iterations-1){
132       outname="RealignObj";
133       outname+=(iter+1);
134       outname.Append(".root");
135       AliAlTrack->WriteRealignObjArray(outname.Data(),AliGeomManager::kSPD1,AliGeomManager::kSSD2);
136     }
137   }
138   
139   if(lastVolid==0){printf("No modules could be realigned \n");return;}
140   printf("End of selecting modules cycle: %d modules selected \n",lastVolid);
141   for(Int_t k=0;k<volIDs2.GetSize();k++){
142     if(volIDs2.At(k)==0)break;
143     layerNum=AliGeomManager::VolUIDToLayer(volIDs2.At(k),modNum);
144     ntVolumeAlign->Fill(layerNum,modNum,volIDs2.At(k));
145   }
146   TFile *f=new TFile("RealignVolNt.root","RECREATE");
147   f->cd();
148   ntVolumeAlign->Write();
149   f->Close();
150
151   timer->Stop();
152   timer->Print();
153   return;
154 }
155
156