PWGUD/multVScentPbPb/AliITSMultRecBg.cxx

   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 ///////////////////////////////////////////////////////////////////////////////
  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
  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(20*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
  52 //_________________________________________________________________
  53 AliITSMultRecBg::AliITSMultRecBg(const AliITSMultRecBg &src)
  54   : AliITSMultReconstructor(src),
  55   fRecType(kData),
  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(),
  66   fInjBuffer(20*kTrNPar)
  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
  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");
 132   if (fRecType==kBgMix && !treeMix) AliFatal("Mixed Mode requested but 2nd RP Tree is missing");
 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();
 544   if (nCurrTr >= fInjBuffer.GetSize()/kTrNPar) fInjBuffer.Set((20+nCurrTr*2)*kTrNPar);
 545   //
 546   Float_t *tracklet = fInjBuffer.GetArray() + nCurrTr*kTrNPar;
 547   Float_t *clPar1=0,*clPar2=0;
 548   //  AliInfo(Form("Size: %d NCurrTr: %d El: %d Winner: %d trackler: %p",fInjBuffer.GetSize(), nCurrTr, nCurrTr*kTrNPar, winner,tracklet));
 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