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