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