ITS UPGRADE
[u/mrichter/AliRoot.git] / ITS / UPGRADE / ITSUpgradeRec / AliITSURecoDet.cxx
1 #include <TGeoVolume.h>
2 #include <TGeoTube.h>
3 #include <TGeoManager.h>
4 #include <TTree.h>
5 #include "AliITSURecoDet.h"
6 #include "AliITSUGeomTGeo.h"
7 #include "AliITSsegmentation.h"
8 #include "AliITSUClusterPix.h"
9 #include "AliITSUReconstructor.h"
10
11
12
13 ClassImp(AliITSURecoDet)
14
15
16 const Char_t* AliITSURecoDet::fgkBeamPipeVolName = "IP_PIPE";
17
18 //______________________________________________________
19 AliITSURecoDet::AliITSURecoDet()
20 :  fNLayers(0)
21   ,fNLayersActive(0)
22   ,fRMax(-1)
23   ,fRMin(-1)
24   ,fRITSTPCRef(-1)
25   ,fLayers(0)
26   ,fLayersActive(0)
27   ,fGeom(0)
28 {
29   // def. c-tor
30 }
31
32 //______________________________________________________
33 AliITSURecoDet::AliITSURecoDet(AliITSUGeomTGeo* geom, const char* name)
34 :  fNLayers(0)
35   ,fNLayersActive(0)
36   ,fRMax(-1)
37   ,fRMin(-1)
38   ,fRITSTPCRef(-1)
39   ,fLayers(0)
40   ,fLayersActive(0)
41   ,fGeom(geom)
42 {
43   // def. c-tor
44   SetNameTitle(name,name);
45   fLayers.SetOwner(kTRUE);        // layers belong to this array
46   fLayersActive.SetOwner(kFALSE); // this one just points on active layers in fLayers
47   Build();
48 }
49
50 //______________________________________________________
51 AliITSURecoDet::~AliITSURecoDet()
52 {
53   // def. d-tor
54   fLayersActive.Clear(); 
55   fLayers.Clear();         // owned!
56 }
57
58 //______________________________________________________
59 void AliITSURecoDet::Print(Option_t* opt) const                       
60 {
61   //print 
62   printf("Detector %s, %d layers, %d active layers\n",GetName(),GetNLayers(),GetNLayersActive());
63   TString opts = opt; opts.ToLower();
64   if (opts.Contains("lr")) {
65     for (int i=0;i<GetNLayers();i++) GetLayer(i)->Print(opt);
66     printf("ITS-TPC matching reference R: %.3f\n",fRITSTPCRef);
67   }
68 }
69
70 //______________________________________________________
71 void AliITSURecoDet::AddLayer(const AliITSURecoLayer* lr)
72 {
73   //add new layer
74   fLayers.AddLast((TObject*)lr);
75   fNLayers++;
76   if (lr->IsActive()) {
77     fLayersActive.AddLast((TObject*)lr);
78     fNLayersActive++;
79   }
80 }
81
82 //______________________________________________________
83 Bool_t AliITSURecoDet::Build()
84 {
85   // build detector from TGeo
86   //
87   if (!fGeom) AliFatal("Geometry interface is not set");
88   int nlr = fGeom->GetNLayers();
89   if (!nlr) AliFatal("No geometry loaded");
90   //
91   // build active ITS layers
92   for (int ilr=0;ilr<nlr;ilr++) {
93     int lrTyp = fGeom->GetLayerChipTypeID(ilr);
94     // name layer according its active id, detector type and segmentation tyoe
95     AliITSURecoLayer* lra = new AliITSURecoLayer(Form("Lr%d%s%d",ilr,fGeom->GetChipTypeName(lrTyp),
96                                                       lrTyp%AliITSUGeomTGeo::kMaxSegmPerChipType),
97                                                  ilr,fGeom);
98     lra->SetPassive(kFALSE);
99     AddLayer(lra);
100   }
101   //
102   // build passive ITS layers
103   //
104   double rMin,rMax,zMin,zMax;
105   // beam pipe
106   TGeoVolume *v = gGeoManager->GetVolume(fgkBeamPipeVolName);
107   AliITSURecoLayer* lrp = 0;
108   if (!v) AliWarning("No beam pipe found in geometry");
109   else {
110     TGeoTube *t=(TGeoTube*)v->GetShape();
111     rMin = t->GetRmin();
112     rMax = t->GetRmax();
113     zMin =-t->GetDz();
114     zMax = t->GetDz();
115     lrp = new AliITSURecoLayer("BeamPipe");
116     lrp->SetRMin(rMin);
117     lrp->SetRMax(rMax);
118     lrp->SetR(0.5*(rMin+rMax));
119     lrp->SetZMin(zMin);
120     lrp->SetZMax(zMax);
121     lrp->SetPassive(kTRUE);
122     AddLayer(lrp);
123     //
124   }
125   //
126   // TPC-ITS wall
127   const AliITSURecoParam* recopar = AliITSUReconstructor::GetRecoParam();
128   if (recopar) {
129     lrp = new AliITSURecoLayer("TPC-ITSwall");
130     lrp->SetRMin(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMin());
131     lrp->SetRMax(AliITSUReconstructor::GetRecoParam()->GetTPCITSWallRMax());
132     lrp->SetR(0.5*(lrp->GetRMin()+lrp->GetRMax()));
133     lrp->SetZMin(-AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
134     lrp->SetZMax( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallZSpanH());
135     lrp->SetMaxStep( AliITSUReconstructor::GetRecoParam()->GetTPCITSWallMaxStep());
136     lrp->SetPassive(kTRUE);
137     AddLayer(lrp);
138   }
139   else {
140     AliWarning("RecoParam is not available, TPC-ITS wall is not set");
141   }
142   //
143   IndexLayers();
144   Print("lr");
145   return kTRUE;
146 }
147
148 //______________________________________________________
149 void AliITSURecoDet::IndexLayers()
150 {
151   // sort and index layers
152   const Double_t kRMargin = 1e-2; // 100 micron margin
153   fLayersActive.Sort();
154   for (int i=0;i<fNLayersActive;i++) GetLayerActive(i)->SetActiveID(i);
155   fLayers.Sort();
156   for (int i=0;i<fNLayers;i++) GetLayer(i)->SetID(i);
157   if (fNLayers>0) {
158     SetRMin(GetLayer(0)->GetRMin()-kRMargin);
159     SetRMax(GetLayer(fNLayers-1)->GetRMax()+kRMargin);
160   }
161   //
162   // define the reference R for ITS/TPC matching: outside of last layer but before TPC materials
163   const double kOffsLastActR = 5.; // offset of reference layer wrt last active R
164   int lastActive = GetLrIDActive(fNLayersActive-1);
165   AliITSURecoLayer* lrA = GetLayer(lastActive); // last active
166   double rref = lrA->GetRMax() + kOffsLastActR;
167   //
168   if (lastActive <  fNLayers-1) { // there are material layers outside ...
169     AliITSURecoLayer* lrL = GetLayer(lastActive+1);
170     if (lrL->GetRMin()<=rref) rref = lrL->GetRMin();
171     if (rref - lrA->GetRMax()<kOffsLastActR) {
172       AliError(Form("The ITS-TPC matching reference R=%.2f is too close to last active R=%.3f",rref,lrA->GetRMax()));
173     }
174   }
175   SetRITSTPCRef(rref);
176   //
177 }
178
179 //______________________________________________________
180 Int_t AliITSURecoDet::FindLastLayerID(Double_t r, int dir) const
181 {
182   // find the last layer which the particle moving in direction dir (1:outward,-1:inward) 
183   // will traverse on its way to radius r 
184   int ilr;
185   //
186   if (dir>0) {
187     for (ilr=0;ilr<fNLayers;ilr++) {
188       AliITSURecoLayer* lr = GetLayer(ilr);
189       if ( r<lr->GetR(-dir) ) break;  // this layer at least entered
190     }
191     return --ilr;  // -1 will correspond to point below the smalles layer
192   }
193   else {
194     for (ilr=fNLayers;ilr--;) {
195       AliITSURecoLayer* lr = GetLayer(ilr);
196       if ( r>lr->GetR(-dir) ) break; // this layer at least entered
197     }
198     ilr++;
199     return ilr<fNLayers ? ilr:-1; // -1 will correspond to point above outer layer
200   }
201   //
202 }
203
204 //______________________________________________________
205 Int_t AliITSURecoDet::FindFirstLayerID(Double_t r, int dir) const
206 {
207   // find the first layer which the particle moving in direction dir (1:outward,-1:inward) 
208   // will traverse starting from radius r 
209   int ilr;
210   //
211   if (dir>0) {
212     for (ilr=0;ilr<fNLayers;ilr++) {
213       AliITSURecoLayer* lr = GetLayer(ilr);
214       if ( r<lr->GetR(dir) ) break;  // this layer at least entered
215     }
216     return ilr<fNLayers ? ilr:-1;  // -1 will correspond to point above outer leayer
217   }
218   else {
219     for (ilr=fNLayers;ilr--;) {
220       AliITSURecoLayer* lr = GetLayer(ilr);
221       if ( r>lr->GetR(dir) ) break; // this layer at least entered
222     }
223     return ilr; // -1 will correspond to point below inner layer
224   }
225   //
226 }
227
228 //______________________________________________________
229 void AliITSURecoDet::CreateClusterArrays()
230 {
231   // create cluster arrays for active layers
232   for (int ilr=0;ilr<fNLayersActive;ilr++) {
233     AliITSURecoLayer*  lr = GetLayerActive(ilr);
234     lr->SetOwnsClusterArray(kTRUE);
235     int tpDet = fGeom->GetLayerChipTypeID(ilr)/AliITSUGeomTGeo::kMaxSegmPerChipType;
236     //
237     if (tpDet == AliITSUGeomTGeo::kChipTypePix) {
238       lr->SetClusters(new TClonesArray(AliITSUClusterPix::Class()));
239     }
240     else {
241       AliFatal(Form("Unknown detector type %d",tpDet));
242     }
243     //
244   }
245   //
246 }
247
248 //_____________________________________________________________________________
249 Int_t AliITSURecoDet::LoadClusters(TTree* treeRP) 
250 {
251   // read clusters from the tree, if it is provided
252   if (!treeRP) return 0;
253   for (int ilr=fNLayersActive;ilr--;) {
254     TBranch* br = treeRP->GetBranch(Form("ITSRecPoints%d",ilr));
255     if (!br) AliFatal(Form("Provided cluster tree does not contain branch for layer %d",ilr));
256     br->SetAddress( GetLayerActive(ilr)->GetClustersAddress() );
257   }
258   return treeRP->GetEntry(0); // we are still in 1 ev/tree mode...
259 }
260
261 Int_t FindIndex(AliITSUClusterPix **parr, Int_t i, AliITSUClusterPix *c) {
262   Int_t b=0;
263   AliITSUClusterPix *cc=parr[b];
264   if (c->Compare(cc) <= 0) return b;
265
266   Int_t e=b+i-1;
267   cc=parr[e];
268   if (c->Compare(cc) >  0) return e+1;
269
270   Int_t m=(b+e)/2;
271   for (; b<e; m=(b+e)/2) {
272     cc=parr[m];
273     if (c->Compare(cc) > 0) b=m+1;
274     else e=m; 
275   }
276   return m;   
277 }
278
279 //______________________________________________________
280 void AliITSURecoDet::SortClusters(AliITSUClusterPix::SortMode_t mode)
281 {
282   // process clsuters according to requested mode
283   AliITSUClusterPix::SetSortMode( mode );
284   for (int ilr=fNLayersActive;ilr--;) {
285       TClonesArray *arr=GetLayerActive(ilr)->GetClusters();
286
287       const Int_t kNcl=arr->GetEntriesFast();
288       AliITSUClusterPix **parr=new AliITSUClusterPix*[kNcl];
289       parr[0]=(AliITSUClusterPix *)arr->UncheckedAt(0);
290       for (Int_t j=1; j<kNcl; j++) {
291           AliITSUClusterPix *c=(AliITSUClusterPix*)arr->UncheckedAt(j);
292           Int_t i=FindIndex(parr,j,c);
293           Int_t k=j-i;
294           memmove(parr+i+1 ,parr+i,k*sizeof(AliITSUClusterPix*));
295           parr[i]=c;
296       }
297
298       TObject **cont=arr->GetObjectRef();
299       for (Int_t i=0; i<kNcl; i++) cont[i]=parr[i];
300       delete[] parr;
301
302   }
303 }
304