]>
Commit | Line | Data |
---|---|---|
a9a39f46 | 1 | #include "AliITSMultRecBg.h" |
2 | #include "AliGeomManager.h" | |
3 | #include "AliMultiplicity.h" | |
4 | #include "../ITS/AliITSgeomTGeo.h" | |
5 | #include <TH2F.h> | |
6 | #include <TTree.h> | |
7 | #include <TRandom.h> | |
8 | #include <TBits.h> | |
9 | #include <TClonesArray.h> | |
10 | ||
c79a7c1c | 11 | /////////////////////////////////////////////////////////////////////////////// |
12 | // // | |
13 | // Class for generating combinatorial backgroung // | |
14 | // for the SPD tracklets // | |
15 | // // | |
16 | // Modes for "insertion", "rotation" and "mixing" are supported // | |
17 | // // | |
18 | /////////////////////////////////////////////////////////////////////////////// | |
19 | ||
20 | ||
a9a39f46 | 21 | ClassImp(AliITSMultRecBg) |
22 | ||
23 | //_________________________________________________________________ | |
24 | AliITSMultRecBg::AliITSMultRecBg() | |
25 | : fRecType(kData), | |
26 | fInjLr(0), | |
27 | fInjStave(0), | |
28 | fInjModule(0), | |
29 | fInjModInStave(0), | |
30 | fInjX(0), | |
31 | fInjZ(0), | |
32 | fInjHitsN(0), | |
33 | fInjScale(1), | |
34 | fInjCurrTrial(0), | |
35 | fInjCluster(), | |
f4c412e5 | 36 | fInjBuffer(20*kTrNPar) |
a9a39f46 | 37 | // |
38 | { | |
39 | // constructor | |
40 | fCreateClustersCopy = kTRUE; | |
41 | for (int i=0;i<2;i++) { | |
42 | fAssociations[i] = 0; | |
43 | fInjSPDOcc[i] = 0; | |
44 | fInjSPDPatt[i] = 0; | |
45 | fInjNTrials[i] = 0; | |
46 | fInjNSuccess[i] = 0; | |
47 | for (int s=4;s--;) fInjModuleClStart[i][s] = fInjModuleClN[i][s] = 0; | |
48 | } | |
49 | // | |
50 | } | |
51 | ||
c79a7c1c | 52 | //_________________________________________________________________ |
ae679e17 | 53 | AliITSMultRecBg::AliITSMultRecBg(const AliITSMultRecBg &src) |
183111ac | 54 | : AliITSMultReconstructor(src), |
ae679e17 | 55 | fRecType(kData), |
c79a7c1c | 56 | fInjLr(0), |
57 | fInjStave(0), | |
58 | fInjModule(0), | |
59 | fInjModInStave(0), | |
60 | fInjX(0), | |
61 | fInjZ(0), | |
62 | fInjHitsN(0), | |
63 | fInjScale(1), | |
64 | fInjCurrTrial(0), | |
65 | fInjCluster(), | |
f4c412e5 | 66 | fInjBuffer(20*kTrNPar) |
c79a7c1c | 67 | // |
68 | { | |
69 | // dummy copy c-tor | |
70 | fCreateClustersCopy = kTRUE; | |
71 | for (int i=0;i<2;i++) { | |
72 | fAssociations[i] = 0; | |
73 | fInjSPDOcc[i] = 0; | |
74 | fInjSPDPatt[i] = 0; | |
75 | fInjNTrials[i] = 0; | |
76 | fInjNSuccess[i] = 0; | |
77 | for (int s=4;s--;) fInjModuleClStart[i][s] = fInjModuleClN[i][s] = 0; | |
78 | } | |
79 | // | |
80 | } | |
81 | ||
a9a39f46 | 82 | //_________________________________________________________________ |
83 | AliITSMultRecBg::~AliITSMultRecBg() | |
84 | { | |
85 | // destructor | |
86 | for (int i=0;i<2;i++) { | |
87 | delete[] fAssociations[i]; | |
88 | delete fInjSPDOcc[i]; | |
89 | delete fInjSPDPatt[i]; | |
90 | for (int s=4;s--;) { | |
91 | delete[] fInjModuleClStart[i][s]; | |
92 | delete[] fInjModuleClN[i][s]; | |
93 | } | |
94 | } | |
95 | // | |
96 | } | |
97 | ||
98 | //____________________________________________________________________ | |
99 | void AliITSMultRecBg::CreateMultiplicityObject() | |
100 | { | |
101 | // create AliMultiplicity object | |
102 | // | |
103 | if (fRecType==kData || fRecType==kBgRot || fRecType==kBgMix) { | |
104 | AliITSMultReconstructor::CreateMultiplicityObject(); | |
105 | } | |
106 | else if (fRecType==kBgInj) { | |
107 | if (fMult) delete fMult; | |
108 | int ntr = GetNInjSuccsess(); | |
109 | TBits fastOrFiredMap(0); | |
110 | fMult = new AliMultiplicity(ntr,0,fNFiredChips[0],fNFiredChips[1],fastOrFiredMap); | |
111 | fMult->SetMultTrackRefs(kTRUE); | |
112 | fMult->SetScaleDThetaBySin2T(fScaleDTBySin2T); | |
113 | fMult->SetDPhiWindow2(fDPhiWindow2); | |
114 | fMult->SetDThetaWindow2(fDThetaWindow2); | |
115 | fMult->SetDPhiShift(fDPhiShift); | |
116 | fMult->SetNStdDev(fNStdDev); | |
117 | // | |
118 | for (int itr=ntr;itr--;) { | |
119 | Float_t *tlInfo = fInjBuffer.GetArray() + itr*kTrNPar; | |
120 | fMult->SetTrackletData(itr,tlInfo); | |
121 | } | |
122 | } | |
123 | // | |
124 | ||
125 | } | |
126 | ||
127 | //_________________________________________________________________________ | |
128 | void AliITSMultRecBg::Run(TTree* tree, Float_t* vtx, TTree* treeMix) | |
129 | { | |
130 | // reconstruct with current settings | |
131 | if (!tree) AliFatal("The RP Tree is missing"); | |
ae679e17 | 132 | if (fRecType==kBgMix && !treeMix) AliFatal("Mixed Mode requested but 2nd RP Tree is missing"); |
a9a39f46 | 133 | if (!vtx) return; |
134 | // | |
135 | if (fRecType==kData) Reconstruct(tree, vtx); | |
136 | else if (fRecType==kBgRot) Reconstruct(tree, vtx); | |
137 | else if (fRecType==kBgMix) ReconstructMix(tree, treeMix, vtx); | |
138 | else if (fRecType==kBgInj) GenInjBgSample(tree,vtx); // if needed, the reco will be done internally | |
139 | else { AliError(Form("Unknown running mode %d",fRecType)); } | |
140 | // | |
141 | CreateMultiplicityObject(); | |
142 | } | |
143 | ||
144 | //_________________________________________________________________________ | |
145 | Bool_t AliITSMultRecBg::PrepareInjBgGenerator(Float_t *vtx) | |
146 | { | |
147 | // prepare histograms/patterns for bg generation | |
148 | // | |
149 | if (!fRecoDone) Reconstruct(fTreeRP,vtx); | |
150 | // | |
151 | float xshift = fSPDSeg.Dx()/2, zshift = fSPDSeg.Dz()/2; | |
152 | int maxRow1 = fSPDSeg.Npx(), maxCol1 = fSPDSeg.Npz(); | |
153 | int nChipsPerLadder = fSPDSeg.GetNumberOfChips(); | |
154 | // | |
155 | fInjCluster.SetLabel(kInjFakeLabel,0); | |
156 | fInjCluster.SetLabel(-2,1); | |
157 | fInjCluster.SetLabel(-2,2); | |
158 | // | |
159 | for (int ilr=0;ilr<2;ilr++) { | |
160 | int nLadPerStave = AliITSgeomTGeo::GetNDetectors(ilr+1); | |
161 | int nStaves = AliITSgeomTGeo::GetNLadders(ilr+1); | |
162 | for (int ild=nLadPerStave;ild--;) { | |
163 | fInjModuleTotClN[ilr][ild] = 0; | |
164 | for (int is=nStaves;is--;) { | |
165 | fInjModuleClStart[ilr][ild][is] = -1; | |
166 | fInjModuleClN[ilr][ild][is] = 0; | |
167 | } | |
168 | } | |
169 | // | |
170 | fInjSPDOcc[ilr]->Reset(); | |
171 | fInjSPDPatt[ilr]->ResetAllBits(); | |
172 | TClonesArray* clusters = fClArr[ilr]; | |
173 | Int_t nclus = fNClustersLay[ilr]; | |
174 | // | |
175 | for (int icl=0;icl<nclus;icl++) { | |
176 | AliITSRecPoint* clus = (AliITSRecPoint*)clusters->UncheckedAt(icl); | |
177 | if (!clus) continue; | |
178 | int ladder = clus->GetDetectorIndex(); // ladder id within the layer | |
179 | int stave = ladder/nLadPerStave; | |
180 | ladder = nLadPerStave-1 - ladder%nLadPerStave; // ladder id within the stave !!!! invert to get human readble histos | |
181 | // the clusters are packed per modules, register the beginning and n.clusters of each module | |
182 | if (fInjModuleClStart[ilr][ladder][stave]<0) fInjModuleClStart[ilr][ladder][stave] = icl; | |
183 | fInjModuleClN[ilr][ladder][stave]++; | |
184 | fInjModuleTotClN[ilr][ladder]++; | |
185 | // | |
186 | int chip = fSPDSeg.GetChipFromLocal(clus->GetDetLocalX(),clus->GetDetLocalZ()); // chip within the ladder | |
187 | if (ilr==1) chip = nChipsPerLadder - 1 - chip; //!!!! invert to get human readble histos | |
188 | chip += ladder*nChipsPerLadder; // chip within the stave | |
189 | ladder %= nLadPerStave; // ladder id within the stave | |
190 | fInjSPDOcc[ilr]->Fill(chip, stave); // register cluster for hit density profile | |
191 | // | |
192 | float xloc = clus->GetDetLocalX()*1e4+xshift; | |
193 | float zloc = clus->GetDetLocalZ()*1e4+zshift; | |
194 | int row,col; | |
195 | fSPDSeg.GetPadIxz(xloc,zloc,row,col); // row,col here stats from 1 | |
196 | row--; | |
197 | col--; | |
198 | // | |
199 | // generate bit pattern according to hit type | |
200 | int npix = clus->GetType(); | |
201 | int nrows = clus->GetNy(); | |
202 | int ncols = clus->GetNz(); | |
203 | float cx=0,cz=0; | |
204 | UInt_t *patt = GenClusterPattern(npix,nrows,ncols,cx,cz); | |
205 | row = int(row - cx); if (row<0) row=0; else if (row+nrows>maxRow1) row = maxRow1-nrows; | |
206 | col = int(col - cz); if (col<0) col=0; else if (col+ncols>maxCol1) col = maxCol1-ncols; | |
207 | for (int icol=ncols;icol--;) { | |
208 | UInt_t hcol = patt[icol]; | |
209 | for (int irow=nrows;irow--;) { | |
210 | if (!(hcol&(1<<irow))) continue; // no hit here | |
211 | int pbit = GetPixelBitL(stave,ladder,col+icol,row+irow); | |
212 | fInjSPDPatt[ilr]->SetBitNumber( pbit ); | |
213 | } | |
214 | } | |
215 | // | |
216 | } // loop over clusters of layer ilr | |
217 | } // loop over layers | |
218 | // | |
219 | if (fNClustersLay[0]==0||fNClustersLay[1]==0) { | |
220 | AliInfo(Form("Trackleting is impossible: Ncl1=%d Ncl2=%d",fNClustersLay[0],fNClustersLay[1])); | |
221 | return kFALSE; | |
222 | } | |
223 | for (int i=0;i<2;i++) { | |
224 | if (fAssociations[i]) delete[] fAssociations[i]; | |
225 | int* tmpi = fAssociations[i] = new Int_t[ fNClustersLay[i] ]; | |
226 | for (int ic=fNClustersLay[i];ic--;) tmpi[ic] = -1; | |
227 | } | |
228 | for (int itr=fNTracklets;itr--;) { | |
229 | Float_t* tracklet = GetTracklet(itr); | |
230 | int cll1 = (int)tracklet[kClID1]; | |
231 | int cll2 = (int)tracklet[kClID2]; | |
232 | fAssociations[0][cll1] = cll2; | |
233 | fAssociations[1][cll2] = cll1; | |
234 | } | |
235 | // | |
236 | fInjBuffer.Set(fNTracklets*kTrNPar); | |
237 | // | |
238 | return kTRUE; | |
239 | // | |
240 | } | |
241 | ||
242 | //_________________________________________________________________________ | |
243 | Bool_t AliITSMultRecBg::GenClusterToInject() | |
244 | { | |
245 | // generate bg cluster on layer lr | |
246 | // | |
247 | int nLadPerStave = AliITSgeomTGeo::GetNDetectors(fInjLr+1); | |
248 | int nChipsPerModule = fSPDSeg.GetNumberOfChips(); | |
249 | // | |
250 | int clID; | |
251 | // | |
252 | //RRR printf("On Lr %d | %d %d\n",fInjLr,fNClustersLay[0],fNClustersLay[1]); | |
253 | do { | |
254 | fInjSPDOcc[fInjLr]->GetRandom2(fInjZ, fInjX); | |
255 | // printf("raw: %f %f\n",fInjZ,fInjX); | |
256 | fInjStave = int(fInjX); | |
257 | int chipInStave = int(fInjZ); // chip in the stave | |
258 | // | |
259 | fInjX = (fInjX - fInjStave)*fSPDSeg.Dx(); // local x coordinate | |
260 | fInjModInStave = chipInStave/nChipsPerModule; | |
261 | fInjModInStave = nLadPerStave - 1 - fInjModInStave; //!!! go from human readable to formal one | |
262 | fInjModule = fInjStave*nLadPerStave + fInjModInStave; // formal module number | |
263 | // | |
264 | fInjZ = (fInjZ-chipInStave)*fSPDSeg.Dz(); // local z coordinate | |
265 | // printf("Z %e X %e | MinSt%d Mod%d\n",fInjZ, fInjX,fInjModInStave,fInjModule); | |
266 | // | |
267 | clID = PickClusterPrototype(fInjLr, fInjModInStave, fInjStave); | |
268 | } while(clID<0); | |
269 | // | |
270 | // printf("clID: %d %d %d %d\n",clID,fNClustersLay[0],fNClustersLay[1],fInjLr); | |
271 | AliITSRecPoint* rClus = (AliITSRecPoint*)fClArr[fInjLr]->UncheckedAt(clID); | |
272 | fInjCluster.SetLayer(fInjLr); | |
273 | fInjCluster.SetType(TMath::Min(kInjMaxNY*kInjMaxNZ,rClus->GetType())); | |
274 | fInjCluster.SetNy(TMath::Min(kInjMaxNY,rClus->GetNy())); | |
275 | fInjCluster.SetNz(TMath::Min(kInjMaxNZ,rClus->GetNz())); | |
276 | fInjCluster.SetDetectorIndex(fInjModule); | |
277 | fInjCluster.SetVolumeId(AliGeomManager::LayerToVolUID(AliGeomManager::kSPD1+fInjLr,fInjModule)); | |
278 | // | |
279 | PlaceInjCluster(); | |
280 | return kTRUE; | |
281 | // | |
282 | } | |
283 | ||
284 | //_________________________________________________________________________ | |
285 | void AliITSMultRecBg::PlaceInjCluster() | |
286 | { | |
287 | // place the injected cluster on the selected module, | |
288 | // avoiding overlaps with real clusters | |
289 | int npix = fInjCluster.GetType(), nrows = fInjCluster.GetNy(), ncols = fInjCluster.GetNz(); | |
290 | Float_t cx=0,cz=0; | |
291 | UInt_t* pattern = GenClusterPattern(npix, nrows, ncols, cx,cz); | |
292 | // | |
293 | // try to embedd on top of real clusters | |
294 | int maxRow1 = fSPDSeg.Npx(), maxCol1 = fSPDSeg.Npz(); | |
295 | int maxRow = maxRow1-1, maxCol = maxCol1-1; | |
296 | TBits &bits = *fInjSPDPatt[fInjLr]; // hits pattern of selected layer | |
297 | int row0=0,col0=0; | |
298 | Bool_t failed; | |
299 | do { | |
300 | failed = kFALSE; | |
301 | fSPDSeg.GetPadIxz(fInjX,fInjZ,row0,col0); | |
302 | row0--; col0--; // row,col here start from 1 | |
303 | if (row0+nrows > maxRow1) row0 = maxRow1-nrows; | |
304 | if (col0+ncols > maxCol1) col0 = maxCol1-ncols; | |
305 | // | |
306 | // printf("Cluster at %f %f: col:%d row:%d\n",fInjZ,fInjX,col0,row0); | |
307 | // | |
308 | // check if generated pattern is mergable with data clusters | |
309 | fInjHitsN = 0; | |
310 | for (int ic=ncols;ic--;) { | |
311 | int colt = col0+ic; | |
312 | for (int ir=nrows;ir--;) { | |
313 | if ( !(pattern[ic]&(1<<ir)) ) continue; // not used | |
314 | int rowt = row0+ir; | |
315 | int bitID = GetPixelBitL(fInjStave,fInjModInStave,colt,rowt); | |
316 | if ( bits.TestBitNumber(bitID) || // is pixel free? | |
317 | (colt>0 && bits.TestBitNumber(bitID-maxRow1)) || // pixel 1 column below | |
318 | (colt<maxCol && bits.TestBitNumber(bitID+maxRow1)) || // pixel 1 column above | |
319 | (rowt>0 && bits.TestBitNumber(bitID-1)) || // pixel 1 row below | |
320 | (rowt<maxRow && bits.TestBitNumber(bitID+1))) // pixel in 1 row above | |
321 | {failed=kTRUE; break;} | |
322 | // ok for this pixel | |
323 | fInjHits[fInjHitsN++] = bitID; | |
324 | } | |
325 | if (failed) break; | |
326 | } | |
327 | if (failed) { // generate new x,z | |
328 | // printf("Conflict found, retry\n"); | |
329 | fInjX = gRandom->Rndm()*fSPDSeg.Dx(); | |
330 | fInjZ = gRandom->Rndm()*fSPDSeg.Dz(); | |
331 | continue; | |
332 | } | |
333 | } while(failed); | |
334 | // | |
335 | // caclulate cluster coordinates | |
336 | float x=0,z=0; | |
337 | fInjX = fInjZ = 0; | |
338 | for (int pix=fInjHitsN;pix--;) { | |
339 | ExtractPixelL(fInjHits[pix], fInjStave, fInjModInStave, col0, row0); | |
340 | fSPDSeg.GetPadCxz(row0+1,col0,x,z); // !!! Note: here row starts from 1 but col from 0!!! | |
341 | fInjX += x; | |
342 | fInjZ += z; | |
343 | } | |
344 | fInjX = (fInjX/fInjHitsN-fSPDSeg.Dx()/2)*1e-4; | |
345 | fInjZ = (fInjZ/fInjHitsN-fSPDSeg.Dz()/2)*1e-4; | |
346 | const TGeoHMatrix *mT2L = AliITSgeomTGeo::GetTracking2LocalMatrix(fInjLr+1,fInjStave+1,fInjModInStave+1); | |
347 | Double_t loc[3]={fInjX,0.,fInjZ},trk[3]={0.,0.,0.}; | |
348 | mT2L->MasterToLocal(loc,trk); | |
349 | // | |
350 | fInjCluster.SetX(0); | |
351 | fInjCluster.SetY(trk[1]); | |
352 | fInjCluster.SetZ(trk[2]); | |
353 | // | |
354 | // printf("ClCoord: %+e %+e %+e\n",fInjCluster.GetX(),fInjCluster.GetY(),fInjCluster.GetZ()); | |
355 | } | |
356 | ||
357 | //_________________________________________________________________________ | |
358 | void AliITSMultRecBg::InitInjBg() | |
359 | { | |
360 | // initialize data for bg injection | |
361 | char buff[100]; | |
362 | for (int ilr=0;ilr<2;ilr++) { | |
363 | sprintf(buff,"occL%d",ilr); | |
364 | int nLaddersStave = AliITSgeomTGeo::GetNDetectors(ilr+1); | |
365 | int nChipsStave = fSPDSeg.GetNumberOfChips()*nLaddersStave; | |
366 | int nStaves = AliITSgeomTGeo::GetNLadders(ilr+1); | |
367 | int nColsStave = fSPDSeg.Npz(); | |
368 | int nRowsStave = fSPDSeg.Npz(); | |
369 | fInjSPDOcc[ilr] = new TH2F(buff,buff,nChipsStave,0,nChipsStave, nStaves,0,nStaves); | |
370 | fInjSPDPatt[ilr] = new TBits(nStaves*nLaddersStave*nColsStave*nRowsStave); | |
371 | for (int is=0;is<nStaves;is++) { | |
372 | for (int il=0;il<nLaddersStave;il++) { | |
373 | fInjModuleClStart[ilr][il] = new Int_t[nStaves]; | |
374 | fInjModuleClN[ilr][il] = new Int_t[nStaves]; | |
375 | } | |
376 | } | |
377 | // | |
378 | } | |
379 | } | |
380 | ||
381 | //_________________________________________________________________________ | |
382 | UInt_t* AliITSMultRecBg::GenClusterPattern(Int_t &npix, Int_t &ny, Int_t &nz, Float_t cy,Float_t &cz) | |
383 | { | |
384 | // generate random pattern for the cluster | |
385 | static UInt_t hitPattern[160]; | |
386 | if (ny>kInjMaxNY) ny = kInjMaxNY; | |
387 | if (nz>kInjMaxNZ) nz = kInjMaxNZ; | |
388 | for (int iz=nz;iz--;) hitPattern[iz]=0; | |
389 | // | |
390 | // handle explicitly easy cases: () is for pixels in the same column ... | |
391 | if (npix==1) hitPattern[0] = 0x1; // type (|) | |
392 | else if (npix==2) { | |
393 | if (nz==1) hitPattern[0] = 0x3; // tpye (||) | |
394 | else hitPattern[0] = hitPattern[1] = 0x1; // tpye (|)(|) | |
395 | } | |
396 | else if (npix==3) { | |
397 | if (nz==1) hitPattern[0] = 0x7; // tpye (|||) | |
398 | else if (ny==1) hitPattern[0] = hitPattern[1] = hitPattern[2] = 0x1; // type (|)(|)(|) | |
399 | else { hitPattern[0] = 0x3; hitPattern[1] = 0x1;} // type (||)(|) | |
400 | } | |
401 | else if (npix==4) { | |
402 | if (nz==1) hitPattern[0] = 0xf; // type (||||) | |
403 | else if (ny==1) hitPattern[0] = hitPattern[1] = hitPattern[2] = hitPattern[3] = 0x1; // type (|)(|)(|)(|) | |
404 | else if (ny==2) { | |
405 | if (nz==2) hitPattern[0] = hitPattern[1] = 0x3; // tpye (||)(||) | |
406 | else { | |
407 | hitPattern[0] = hitPattern[1] = hitPattern[2] = 0x1; // type (||)(|)(|) or (|)(||)(|) | |
408 | hitPattern[gRandom->Rndm()>0.5 ? 0 : 1] |= 0x2; | |
409 | } | |
410 | } | |
411 | else { | |
412 | hitPattern[0] = 0x7; | |
413 | hitPattern[1] = (gRandom->Rndm()>0.8) ? 0x1 : 0x2; // type (|||) + (_|_) or (|||) + (|__) | |
414 | } | |
415 | } | |
416 | // more complex topologies | |
417 | else { | |
418 | UInt_t mask = 0xffffffff>>(32-ny); | |
419 | for (int i=nz;i--;) hitPattern[i] = mask; | |
420 | int toSup = ny*nz - npix; | |
421 | int trial = toSup*5; | |
422 | int colToTouch = nz-1; // at least 1 column should not be touched | |
423 | while(toSup>0 && (trial--)>0) { // suppress random pixel until needed npix, ny and nz is obtained | |
424 | int col = gRandom->Integer(nz); | |
425 | if (hitPattern[col]==0x1) continue; // don't lose this column | |
426 | if (hitPattern[col]==mask) { | |
427 | if (!colToTouch) continue; // this is the only remaining column with ny rows hit | |
428 | colToTouch--; // there are other columns of ny rows hit, may touch it | |
429 | } | |
430 | hitPattern[col] >>= 1; | |
431 | toSup--; | |
432 | } | |
433 | if (toSup) npix += toSup; // failed to suppress exact number of pixels | |
434 | } | |
435 | // | |
436 | cy = cz = 0; // get centroid | |
437 | for (int col=nz;col--;) { | |
438 | int npxr = 0; | |
439 | for (int row=ny;row--;) if (hitPattern[col] & (0x1<<row)) {npxr++; cy+=(1+row);} | |
440 | cz += npxr*(1+col); | |
441 | } | |
442 | cz = cz/npix - 0.5; | |
443 | cy = cy/npix - 0.5; | |
444 | // | |
445 | return (UInt_t*)&hitPattern[0]; | |
446 | } | |
447 | ||
448 | //_________________________________________________________________________ | |
449 | Int_t AliITSMultRecBg::PickClusterPrototype(Int_t lr, Int_t ladInStave, Int_t stave2avoid) | |
450 | { | |
451 | // pick random cluster to use as a prototype for injection. It should come | |
452 | // from ladders at given Z belt and should not be from the stave2avoid | |
453 | // pick random cluster to use as a prototype for injection. It should come | |
454 | // from ladders at given Z belt and should not be from the stave2avoid | |
455 | static int tried = 0; | |
456 | static int ladInStaveOrig = 0; | |
457 | // | |
458 | int ncl = fInjModuleTotClN[lr][ladInStave]; | |
459 | if (stave2avoid>=0) ncl -= fInjModuleClN[lr][ladInStave][stave2avoid]; | |
460 | if (ncl<1) { | |
461 | int totLad = AliITSgeomTGeo::GetNDetectors(lr+1); | |
462 | if (!tried) ladInStaveOrig = ladInStave; // before starting the resque solution, save the orig.ladder | |
463 | if (++tried >= totLad) { // failed to find cluster avoiding a stave2avoid | |
464 | tried = 0; | |
465 | return PickClusterPrototype(lr,ladInStaveOrig,-1); | |
466 | } | |
467 | int useLad = ladInStave+1; | |
468 | if (useLad>=AliITSgeomTGeo::GetNDetectors(lr+1)) useLad = 0; | |
469 | return PickClusterPrototype(lr,useLad,stave2avoid); // look in the neigbouring ladder | |
470 | } | |
471 | // | |
472 | int pick = gRandom->Integer(ncl); | |
473 | int nst = AliITSgeomTGeo::GetNLadders(lr+1); | |
474 | int stave = 0; | |
475 | for (stave=0;stave<nst;stave++) { | |
476 | if (stave==stave2avoid) continue; | |
477 | if (pick<fInjModuleClN[lr][ladInStave][stave]) break; | |
478 | pick -= fInjModuleClN[lr][ladInStave][stave]; | |
479 | } | |
480 | // | |
481 | tried = 0; | |
482 | return fInjModuleClStart[lr][ladInStave][stave]+pick; | |
483 | } | |
484 | ||
485 | //_________________________________________________________________________ | |
486 | Int_t AliITSMultRecBg::SearchInjTracklet(const Float_t *vtx) | |
487 | { | |
488 | // reconstruct tracklets which may involve injected cluster | |
489 | // fake arrays to be used for injected cluster in MultReco | |
490 | Float_t clustersLayInj[kClNPar]; | |
491 | Int_t detectorIndexClustersLayInj[1]; | |
492 | Bool_t overlapFlagClustersLayInj[1]; | |
493 | UInt_t usedClusLayInj[1]; | |
494 | // | |
495 | Bool_t kUseOrig = kFALSE;//kTRUE; | |
496 | if (kUseOrig) { | |
497 | // to try with unused clusters | |
498 | if ( fAssociations[fInjLr][fInjCurrTrial]>=0 ) return 0; // associated | |
499 | float *origCl = GetClusterOfLayer(fInjLr,fInjCurrTrial); | |
500 | for (int i=0;i<kClNPar;i++) clustersLayInj[i] = origCl[i]; | |
501 | } | |
502 | else { | |
503 | // >> fill cluster data: equivavlent of data fetched in AliITSMultReconstructor::LoadClusterArrays | |
504 | detectorIndexClustersLayInj[0] = fInjCluster.GetDetectorIndex(); | |
505 | overlapFlagClustersLayInj[0] = kFALSE; | |
506 | usedClusLayInj[0] = 0; | |
507 | fInjCluster.GetGlobalXYZ( clustersLayInj ); | |
508 | for (int i=3;i--;) clustersLayInj[kClMC0+i] = fInjCluster.GetLabel(i); | |
509 | ClusterPos2Angles(clustersLayInj, vtx); // convert to angles | |
510 | } | |
511 | // | |
512 | // compare injected cluster with real ones | |
513 | int partnerLr = 1 - fInjLr; | |
514 | ||
515 | double bestChiUsed = fNStdDev*2, bestChiFree = fNStdDev*2; | |
516 | int bestPartnerUsed=-1,bestPartnerFree=-1; | |
517 | // | |
518 | for (int icl=fNClustersLay[partnerLr];icl--;) { | |
519 | Float_t *partnerCl = GetClusterOfLayer(partnerLr,icl); | |
520 | Double_t dTheta = clustersLayInj[kClTh] - partnerCl[kClTh]; | |
521 | Double_t dPhi = clustersLayInj[kClPh] - partnerCl[kClPh]; | |
522 | if (dPhi>TMath::Pi()) dPhi=2.*TMath::Pi()-dPhi; // take into account boundary condition | |
523 | Float_t d = CalcDist(dPhi,dTheta, fInjLr==0 ? clustersLayInj[kClTh]:partnerCl[kClTh]); | |
524 | if (d>fNStdDev) continue; | |
525 | // | |
526 | int competitor = fAssociations[partnerLr][icl]; // is the cluster of partner layer already used by some tracklet? | |
527 | if (competitor>=0) { if (d<bestChiUsed && d < fMinDists[ partnerLr==1 ? icl : competitor ]) { bestChiUsed = d; bestPartnerUsed = icl; } } | |
528 | else { if (d<bestChiFree) { bestChiFree = d; bestPartnerFree = icl; } } | |
529 | } | |
530 | // | |
531 | int winner = -1; | |
532 | // | |
533 | if (bestPartnerUsed>=0) { | |
534 | if (bestPartnerFree>=0) { winner = bestChiFree<bestChiUsed ? bestPartnerFree : bestPartnerUsed;} // shall we subtract old real tracklet if the winner cluster is used one? | |
535 | else { winner = bestPartnerUsed; } | |
536 | } | |
537 | else winner = bestPartnerFree; | |
538 | ||
539 | // printf("%d\n",winner); | |
540 | ||
541 | if (winner<0) return 0; | |
542 | // | |
543 | int nCurrTr = GetNInjSuccsess(); | |
f4c412e5 | 544 | if (nCurrTr >= fInjBuffer.GetSize()/kTrNPar) fInjBuffer.Set((20+nCurrTr*2)*kTrNPar); |
a9a39f46 | 545 | // |
546 | Float_t *tracklet = fInjBuffer.GetArray() + nCurrTr*kTrNPar; | |
547 | Float_t *clPar1=0,*clPar2=0; | |
f4c412e5 | 548 | // AliInfo(Form("Size: %d NCurrTr: %d El: %d Winner: %d trackler: %p",fInjBuffer.GetSize(), nCurrTr, nCurrTr*kTrNPar, winner,tracklet)); |
a9a39f46 | 549 | if (fInjLr) { |
550 | clPar1 = GetClusterOfLayer(partnerLr,winner); | |
551 | clPar2 = clustersLayInj; | |
552 | tracklet[kClID1] = winner; | |
553 | tracklet[kClID2] = -1; | |
554 | } | |
555 | else { | |
556 | clPar2 = GetClusterOfLayer(partnerLr,winner); | |
557 | clPar1 = clustersLayInj; | |
558 | tracklet[kClID1] = -1; | |
559 | tracklet[kClID2] = winner; | |
560 | } | |
561 | tracklet[kTrTheta] = clPar1[kClTh]; // use the theta from the clusters in the first layer | |
562 | tracklet[kTrPhi] = clPar1[kClPh]; // use the phi from the clusters in the first layer | |
563 | tracklet[kTrDPhi] = clPar1[kClPh] - clPar2[kClPh]; // store the difference between phi1 and phi2 | |
564 | // | |
565 | // define dphi in the range [0,pi] with proper sign (track charge correlated) | |
566 | if (tracklet[kTrDPhi] > TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]-2.*TMath::Pi(); | |
567 | if (tracklet[kTrDPhi] < -TMath::Pi()) tracklet[kTrDPhi] = tracklet[kTrDPhi]+2.*TMath::Pi(); | |
568 | tracklet[kTrDTheta] = clPar1[kClTh] - clPar2[kClTh]; // store the theta1-theta2 | |
569 | tracklet[kTrLab1] = clPar1[kClMC0]; | |
570 | tracklet[kTrLab2] = clPar2[kClMC0]; | |
571 | // | |
572 | // printf("Got Tracklet from lr%d: %f %f\n",fInjLr,fInjNSuccess[0],fInjNSuccess[1]); | |
573 | // | |
574 | return 1; | |
575 | // | |
576 | } | |
577 | ||
578 | //_________________________________________________________________________ | |
579 | void AliITSMultRecBg::GenInjBgSample(TTree* treeRP, Float_t *vtx) | |
580 | { | |
581 | // generate a sample of tracklets from injected bg | |
582 | // | |
583 | SetTreeRP(treeRP); | |
584 | InitInjBg(); | |
585 | if (!PrepareInjBgGenerator(vtx)) return; | |
586 | // | |
587 | fInjNSuccess[0] = fInjNSuccess[1] = 0; | |
588 | for (int i=0;i<2;i++) { | |
589 | fInjLr = i; | |
590 | fInjNTrials[i] = TMath::Max(1,int(fInjScale*fNClustersLay[i])); | |
591 | for (fInjCurrTrial=0;fInjCurrTrial<fInjNTrials[i];fInjCurrTrial++) { | |
592 | if (!GenClusterToInject()) break; | |
593 | fInjNSuccess[i] += SearchInjTracklet(vtx); | |
594 | } | |
595 | } | |
596 | printf("Successes/Trials: %d/%d %d/%d\n",fInjNSuccess[0],fInjNTrials[0],fInjNSuccess[1],fInjNTrials[1]); | |
597 | // | |
598 | } | |
599 |