]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/multVScentPbPb/AliITSMultRecBg.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGUD / multVScentPbPb / AliITSMultRecBg.cxx
CommitLineData
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 21ClassImp(AliITSMultRecBg)
22
23//_________________________________________________________________
24AliITSMultRecBg::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 53AliITSMultRecBg::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//_________________________________________________________________
83AliITSMultRecBg::~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//____________________________________________________________________
99void 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//_________________________________________________________________________
128void 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//_________________________________________________________________________
145Bool_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//_________________________________________________________________________
243Bool_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//_________________________________________________________________________
285void 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//_________________________________________________________________________
358void 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//_________________________________________________________________________
382UInt_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//_________________________________________________________________________
449Int_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//_________________________________________________________________________
486Int_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
032b4329 523 Float_t d = CalcDist(dPhi,fInjLr==0 ? -dTheta:dTheta, fInjLr==0 ? clustersLayInj[kClTh]:partnerCl[kClTh]);
a9a39f46 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//_________________________________________________________________________
579void 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