EVGEN/AliGenHBTosl.cxx

   1 #include "AliGenHBTosl.h"
   2 //__________________________________________________________
   3 /////////////////////////////////////////////////////////////
   4 //                                                         //
   5 //  class AliGenHBTosl                                     //
   6 //                                                         //
   7 //  Genarator simulating particle correlations             //
   8 //                                                         //
   9 //  The main idea of the generator is to produce particles //
  10 //  according to some distribution of two particle         //
  11 //  property. In HBT they are qout,qsie and qlong.         //
  12 //  In order to be able to generate signal that produces   //
  13 //  given two particle correlation background must be      //
  14 //  known before in order to produce the shape of signal   //
  15 //  to randomize given distribution from.                  //
  16 //                                                         //
  17 //  The generator works as follows:                        //
  18 //  1. Coarse Background (fQCoarseBackground) is generated //
  19 //     ade  from the particles                             //
  20 //     given by the external generator (variable           //
  21 //     fGenerator) by the mixing technique.                //
  22 //  2. Coarse signal is prduced by multiplying Coarse      //
  23 //     background by a required function                   //
  24 //     See method FillCoarseSignal                         //
  25 //  3. Signal is randomized out of the coarse signal       //
  26 //     histogram (two particle property). First particle   //
  27 //     is taken from the external generator, and the       //
  28 //     second one is CALCULATED on the basis of the first  //
  29 //     one and the two particle property (qout,qside,qlong)//
  30 //     Background is made by the mixing out of the         //
  31 //     genereted signal events.                            //
  32 //     This step is cotinued up to the moment signal       //
  33 //     histogram has enough statistics (data member        //
  34 //     fMinFill)                                           //
  35 //     See method StartSignalPass1()                       //
  36 //  4. chi is calculated for each bin (chiarray variqable) //
  37 //     (not the chi2 because sign is important)            //
  38 //     Two particle prioperty                              //
  39 //     (qout,qside,qlong) is chosen at the points that     //
  40 //     chi is the smallest. First particle is taken from   //
  41 //     the the external generator (fGenerator) and second's /
  42 //     momenta are caclulated out of their momenta and     //
  43 //     (qout,qside,qlong). Background is updated           //
  44 //     continuesely for all the events. This step is       //
  45 //     continued until stability conditions are fullfiled  //
  46 //     or maximum number of iteration is reached.          //
  47 //  5. The same as step 4 but events are stored.           //
  48 //                                                         //
  49 ////////////////////////////////////////////////////////////
  50
  51 #include <TCanvas.h>
  52
  53
  54 #include <TH3D.h>
  55 #include <TList.h>
  56 #include <TPDGCode.h>
  57 #include <TParticle.h>
  58 #include <AliStack.h>
  59 #include <TMath.h>
  60 #include <TVector3.h>
  61 #include <TStopwatch.h>
  62 #include <TFile.h>
  63
  64 #include "AliGenCocktailAfterBurner.h"
  65 #include "AliGeVSimParticle.h"
  66 #include "AliGenGeVSim.h"
  67 #include "AliGenHIJINGpara.h"
  68
  69
  70 /***********************************************************/
  71 ClassImp(AliGenHBTosl)
  72
  73 AliGenHBTosl::AliGenHBTosl():
  74  AliGenerator(),
  75  fQCoarseBackground(0x0),
  76  fQCoarseSignal(0x0),
  77  fQSignal(0x0),
  78  fQBackground(0x0),
  79  fQSecondSignal(0x0),
  80  fQSecondBackground(0x0),
  81  fQRange(0.06),
  82  fQNBins(60),
  83  fGenerator(0x0),
  84  fBufferSize(5),
  85  fNBinsToScale(Int_t(fQNBins*0.1)),
  86  fDebug(0),
  87  fMaxIterations(1000),
  88  fMaxChiSquereChange(0.01),
  89  fMaxChiSquerePerNDF(1.5),
  90  fQRadius(8.0),
  91  fPID(kPiPlus),
  92  fSamplePhiMin(-0.01),
  93  fSamplePhiMax(TMath::TwoPi()+0.01),
  94  fSignalRegion(0.0),
  95  fMinFill(1000),
  96  fLogFile(0x0)
  97 {
  98 //default constructor
  99 }
 100 /***********************************************************/
 101
 102 AliGenHBTosl::AliGenHBTosl(Int_t n, Int_t pid):
 103  AliGenerator(n),
 104  fQCoarseBackground(0x0),
 105  fQCoarseSignal(0x0),
 106  fQSignal(0x0),
 107  fQBackground(0x0),
 108  fQSecondSignal(0x0),
 109  fQSecondBackground(0x0),
 110  fQRange(0.06),
 111  fQNBins(60),
 112  fGenerator(0x0),
 113  fBufferSize(5),
 114  fNBinsToScale(Int_t(fQNBins*0.1)),
 115  fDebug(0),
 116  fSignalShapeCreated(kFALSE),
 117  fMaxIterations(1000),
 118  fMaxChiSquereChange(0.01),
 119  fMaxChiSquerePerNDF(1.5),
 120  fQRadius(8.0),
 121  fPID(pid),
 122  fSamplePhiMin(-0.01),
 123  fSamplePhiMax(TMath::TwoPi()+0.01),
 124  fSignalRegion(0.0),
 125  fMinFill(1000),
 126  fLogFile(0x0)
 127 {
 128 //default constructor
 129 }
 130 /***********************************************************/
 131
 132 AliGenHBTosl::~AliGenHBTosl()
 133 {
 134 //destructor
 135  delete fQCoarseSignal;
 136  delete fQCoarseBackground;
 137  delete fQSignal;
 138  delete fQBackground;
 139  delete fGenerator;
 140  delete fQSecondSignal;
 141  delete fQSecondBackground;
 142  delete fLogFile;
 143 }
 144 /***********************************************************/
 145
 146 void AliGenHBTosl::Init()
 147 {
 148   //Initializes generator
 149   if (fGenerator == 0x0)
 150    {
 151
 152      AliGenHIJINGpara* bkggen = new AliGenHIJINGpara(fNpart*4);
 153      fGenerator = bkggen;
 154
 155 /*
 156      AliGenGeVSim * gevsim = new AliGenGeVSim(0.0);
 157      AliGeVSimParticle* kplus = new AliGeVSimParticle(fPID,1,fNpart, 0.17, 0.9);
 158      gevsim->AddParticleType(kplus);
 159
 160      fGenerator = gevsim;
 161 */
 162
 163 /*
 164
 165     AliMevSimConfig *c = new AliMevSimConfig(1);
 166     c->SetRectPlane(1);                                 // reaction plane control, model 4
 167     c->SetGrid(80,80);
 168
 169     AliGenMevSim *mevsim = new AliGenMevSim(c);
 170     mevsim->SetPtRange(0.001, 3);
 171     mevsim->SetMomentumRange(0.1, 3);
 172     mevsim->SetTrackingFlag(0);
 173     mevsim->SetOrigin(0.0, 0.0, 0.0);
 174     mevsim->SetSigma(0.0, 0.0, 0.0);
 175     AliMevSimParticle *kplus = new AliMevSimParticle(kKPlus, fNpart, 0, 0.25, 0.0, 2, 0.15, 0.0, 0.0 );
 176     mevsim->AddParticleType(kplus);
 177     fGenerator = mevsim;
 178 */
 179
 180     fGenerator->SetOrigin(fOrigin[0],fOrigin[1],fOrigin[2]);
 181     static const Double_t kDegToRadCF = 180./TMath::Pi();
 182     fGenerator->SetMomentumRange(fPtMin,fPtMax);
 183     fGenerator->SetPhiRange(kDegToRadCF*fPhiMin,kDegToRadCF*fPhiMax);
 184     fGenerator->SetYRange(fYMin,fYMax);
 185     fGenerator->SetThetaRange(kDegToRadCF*fThetaMin,kDegToRadCF*fThetaMax);
 186     fGenerator->Init();
 187
 188    }
 189
 190 //  fQCoarseBackground =  new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 191 //  fQCoarseSignal =  new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 192 //  fQCoarseBackground->Sumw2();
 193 //  fQCoarseSignal->Sumw2();
 194
 195   fQSignal =  new TH3D("fQSignal1","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 196   fQBackground =  new TH3D("fQBackground1","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 197
 198   fQSecondSignal =  new TH3D("fQSignal2","fQSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 199   fQSecondBackground =  new TH3D("fQBackground2","fQBackground",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 200
 201   fQSignal->Sumw2();
 202   fQBackground->Sumw2();
 203   fQSecondSignal->Sumw2();
 204   fQSecondBackground->Sumw2();
 205
 206   fLogFile = new ofstream("BadEvent",ios::out);
 207
 208 }
 209 /***********************************************************/
 210
 211 void AliGenHBTosl::Generate()
 212 {
 213  //the main method
 214
 215  ofstream& logfile = *fLogFile;
 216  logfile<<"Generate"<<"Attempts to generate "<<fNpart<<" particles.";
 217
 218
 219  if (fStackBuffer == 0x0) fStackBuffer = new TList();
 220  //Here is initialization level
 221  if (fSignalShapeCreated == kFALSE)
 222   {
 223    TH3D *hs = 0x0, *hb = 0x0;
 224    TFile* file;
 225
 226    file = TFile::Open("QTSignal.root");
 227    if (file)
 228     {
 229      hs = (TH3D*)file->Get("fQSignal1");
 230      if (hs) hs->SetDirectory(0x0);
 231     }
 232    delete file;
 233
 234    file = TFile::Open("QTBackground.root");
 235    if (file)
 236     {
 237      hb = (TH3D*)file->Get("fQBackground1");
 238      if (hb) hb->SetDirectory(0x0);
 239     }
 240    delete file;
 241
 242    if (hs && hb)
 243     {
 244       Info("Generate","**********************************");
 245       Info("Generate","Found starting histograms in files");
 246       Info("Generate","**********************************");
 247       delete fQSignal;
 248       delete fQBackground;
 249       fQSignal = hs;
 250       fQBackground = hb;
 251     }
 252    else
 253     {
 254       TH3D *cs = 0x0, *cb = 0x0;
 255       file = TFile::Open("QTCoarseBackground.root");
 256       if (file)
 257        {
 258         cb = (TH3D*)file->Get("fQCoarseBackground");
 259         if (cb) cb->SetDirectory(0x0);
 260        }
 261       delete file;
 262
 263       file = TFile::Open("QTCoarseSignal.root");
 264       if (file)
 265        {
 266         cs = (TH3D*)file->Get("fQCoarseSignal");
 267         if (cs) cs->SetDirectory(0x0);
 268        }
 269       delete file;
 270
 271       if (cs && cb)
 272        {
 273
 274          Info("Generate","Got Coarse signal and bkg from files");
 275          delete fQCoarseBackground;
 276          delete fQCoarseSignal;
 277          fQCoarseSignal = cs;
 278          fQCoarseBackground = cb;
 279        }
 280       else
 281        {
 282          if (cb)
 283            {
 284              Info("Generate","Got Coarse bkg from file");
 285              delete fQCoarseBackground;
 286              fQCoarseBackground = cb;
 287            }
 288          else
 289            {
 290              fQCoarseBackground =  new TH3D("fQCoarseBackground","",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 291              fQCoarseBackground->Sumw2();
 292
 293              FillCoarse();      //create coarse background - just to know the spectrum
 294            }
 295
 296          fQCoarseSignal =  new TH3D("fQCoarseSignal","fQCoarseSignal",fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange, fQNBins,-fQRange,fQRange);
 297          fQCoarseSignal->Sumw2();
 298          FillCoarseSignal();//create first coarse signal by brutal multplication coarse background and required function shape
 299        }
 300
 301       StartSignal();     //Initilizes the stack that is used for generation
 302     }
 303    fSignalShapeCreated = kTRUE;
 304   }
 305
 306  AliStack* stack = RotateStack();
 307
 308  AliStack* genstack = fGenerator->GetStack();
 309  if (genstack == 0x0)
 310   {
 311    genstack = new AliStack(fNpart);
 312    fGenerator->SetStack(genstack);
 313   }
 314  else
 315   {
 316    genstack->Reset();
 317   }
 318
 319  fGenerator->Generate();
 320  Int_t j = 0, ntr = 0;
 321  if ( genstack->GetNtrack() < fNpart/2)
 322   {
 323     Warning("Generate","************************************************************");
 324     Warning("Generate","Generator generated (%d) less particles then expected (%d).",
 325              stack->GetNtrack(),fNpart/2);
 326     Warning("Generate","************************************************************");
 327   }
 328
 329  TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
 330  work->SetDirectory(0x0);
 331  work->Sumw2();
 332
 333  Double_t*** chiarray = new Double_t** [fQNBins+1];
 334  Double_t*** sigarray = new Double_t** [fQNBins+1];
 335
 336  for (Int_t i = 1; i<=fQNBins; i++)
 337    {
 338       chiarray[i] =  new Double_t* [fQNBins+1];
 339       sigarray[i] =  new Double_t* [fQNBins+1];
 340
 341       for (Int_t k = 1; k<=fQNBins; k++)
 342        {
 343          chiarray[i][k] =  new Double_t [fQNBins+1];
 344          sigarray[i][k] =  new Double_t [fQNBins+1];
 345        }
 346    }
 347
 348
 349  Double_t scale = Scale(fQSignal,fQBackground);
 350  work->Divide(fQSignal,fQBackground,scale);
 351
 352  Double_t binwdh = work->GetBinWidth(1)/2.;
 353
 354  for (Int_t k = 1; k<=fQNBins; k++)
 355    {
 356     Double_t z = work->GetZaxis()->GetBinCenter(k);
 357     for (Int_t j = 1; j<=fQNBins; j++)
 358       {
 359        Double_t y = work->GetYaxis()->GetBinCenter(j);
 360        for (Int_t i = 1; i<=fQNBins; i++)
 361          {
 362               sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
 363               Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qinv)
 364               Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
 365               Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
 366               chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
 367          }
 368        }
 369     }
 370
 371  char msg[1000];
 372  logfile<<endl;
 373  sprintf(msg,"\n");
 374  Int_t middlebin = fQNBins/2;
 375
 376  for (Int_t k = middlebin-5; k < middlebin+5; k++)
 377    {
 378      Double_t tx = work->GetXaxis()->GetBinCenter(30);
 379      Double_t ty = work->GetYaxis()->GetBinCenter(30);
 380      Double_t tz = work->GetZaxis()->GetBinCenter(k);
 381      sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
 382      logfile<<msg;
 383    }
 384  logfile<<endl;
 385
 386  for (Int_t k = middlebin-5; k < middlebin+5; k++)
 387   {
 388     sprintf(msg,"% 6.5f ",work->GetBinContent(30,30,k));
 389     logfile<<msg;
 390   }
 391  logfile<<endl;
 392
 393  for (Int_t k = middlebin-5; k < middlebin+5; k++)
 394   {
 395     sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
 396     logfile<<msg;
 397   }
 398  logfile<<endl;
 399
 400  TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
 401  TParticle* second = &particle;
 402
 403  Bool_t shortloop = kTRUE;
 404  Int_t sc = 0;//security check against infinite loop
 405
 406  while ( (ntr+1) < fNpart)
 407    {
 408      Int_t xmax = 1;
 409      Int_t ymax = 1;
 410      Int_t zmax = 1;
 411      Double_t qout;
 412      Double_t qside;
 413      Double_t qlong;
 414
 415
 416      Int_t loopmin;
 417      Int_t loopmax;
 418
 419      if (shortloop)
 420       {
 421         shortloop = kFALSE;
 422         loopmax = fQNBins;
 423         loopmin = 1;
 424       }
 425      else
 426       {
 427         shortloop = kTRUE;
 428         loopmax = fQNBins/2+fQNBins/4;
 429         loopmin = fQNBins/2-fQNBins/4;
 430       }
 431
 432
 433      for (Int_t k = loopmin; k <=loopmax; k++ )
 434       {
 435        qlong = work->GetZaxis()->GetBinCenter(k);
 436        for (Int_t j = loopmin; j<=loopmax; j++)
 437          {
 438           qside = work->GetYaxis()->GetBinCenter(j);
 439           for (Int_t i = loopmin; i<=loopmax; i++)
 440             {
 441               qout  = work->GetXaxis()->GetBinCenter(i);
 442               if (chiarray[xmax][ymax][zmax] < chiarray[i][j][k])
 443                {
 444                  xmax = i;
 445                  ymax = j;
 446                  zmax = k;
 447                }
 448
 449 //              Double_t qdist = TMath::Sqrt(qout*qout + qside*qside + qlong*qlong);
 450
 451 //              Double_t fact = chiarray[i][j][k];//chiarray is chi2
 452 //              if (fact > work->GetBinError(i,j,k))//if differece between what we want and
 453 //               {                                  //what we have is bigger than stat. error
 454 //                 xmax = i;                        //we force to fill that bin
 455 //                 ymax = j;
 456 //                 zmax = k;
 457 //                 break;
 458 //               }
 459             }
 460          }
 461       }
 462      Double_t qlongc = work->GetZaxis()->GetBinCenter(zmax);
 463      Double_t qsidec = work->GetYaxis()->GetBinCenter(ymax);
 464      Double_t qoutc  = work->GetXaxis()->GetBinCenter(xmax);
 465
 466
 467      sprintf(msg,"Generate Fill bin chi2(%d,%d,%d)=%f",xmax,ymax,zmax,chiarray[xmax][ymax][zmax]);
 468      logfile<<msg;
 469
 470      qout  = gRandom->Uniform(qoutc -binwdh, qoutc +binwdh);
 471      qside = gRandom->Uniform(qsidec-binwdh, qsidec+binwdh);
 472      qlong = gRandom->Uniform(qlongc-binwdh, qlongc+binwdh);
 473
 474      TParticle* first = 0;
 475      while (j < genstack->GetNtrack())
 476       {
 477         TParticle* tmpp = genstack->Particle(j++);
 478         if (tmpp->GetPdgCode() == fPID)
 479          {
 480            if (CheckParticle(tmpp,0x0,stack) == kFALSE)
 481             {
 482               first = tmpp;
 483               break;
 484             }
 485          }
 486       }
 487
 488      if (first == 0x0)
 489       {
 490         if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
 491         break;
 492       }
 493
 494      //Here put the check if this particle do not fall into signal region with other partticle
 495
 496       Int_t retval = GetThreeD(first,second,qout,qside,qlong);
 497       if (retval)
 498        {
 499          //Info("StartSignal","Can not find momenta for this OSL and particle");
 500          continue;
 501        }
 502     //in case this particle is falling into signal area with another
 503     //particle we take a another pair
 504     //it can intruduce artificial correlations
 505      Bool_t checkresult = CheckParticle(second,first,stack);
 506      if ( checkresult  && (sc < 10) )
 507       {
 508         sc++;
 509         continue;
 510       }
 511      sc = 0;
 512
 513      //Put on output stack
 514      SetTrack(first,ntr);
 515      SetTrack(second,ntr);
 516
 517      //Put on internal stack
 518      Int_t etmp;
 519      SetTrack(first,etmp,stack);
 520      SetTrack(second,etmp,stack);
 521
 522      Double_t y = GetQOutQSideQLongCorrTheorValue(qoutc,qsidec,qlongc);
 523
 524      sigarray[xmax][ymax][zmax] ++;
 525      chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
 526      chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
 527
 528     }
 529
 530  Mix(fStackBuffer,fQBackground,fQSecondSignal); //upgrate background
 531  Mix(stack,fQSignal,fQSecondBackground); //upgrate signal
 532
 533  delete work;
 534
 535  for (Int_t i = 1; i<=fQNBins; i++)
 536    {
 537      for (Int_t k = 1; k<=fQNBins; k++)
 538       {
 539         delete [] chiarray[i][k];
 540         delete [] sigarray[i][k];
 541       }
 542      delete [] chiarray[i];
 543      delete [] sigarray[i];
 544    }
 545  delete [] chiarray;
 546  delete [] sigarray;
 547 }
 548 /***********************************************************/
 549
 550 void AliGenHBTosl::GetOneD(TParticle* first, TParticle* second,Double_t qinv)
 551 {
 552 //deprecated method that caclulates momenta of the second particle
 553 // out of qinv and the first particle
 554 //first particle is rotated that only X is non-zero
 555
 556
 557   Double_t m = first->GetMass();
 558   Double_t msqrd = m*m;
 559   Double_t fourmassSquered = 4.*msqrd;
 560
 561   //Condition that R must fullfill to be possible to have qinv less smaller then randomized
 562 //  Double_t rRange = qinv*TMath::Sqrt(qinv*qinv + fourmassSquered)/fourmassSquered;
 563 //  Double_t r = gRandom->Uniform(rRange);
 564
 565   Double_t r = gRandom->Uniform(qinv);
 566   Double_t phi = gRandom->Uniform(TMath::TwoPi());
 567
 568   Double_t firstPx = first->P();//first particle is rotated that only X is non-zero  thus P==Px
 569   Double_t px = 2.*msqrd*firstPx + firstPx*qinv*qinv;
 570   Double_t temp = qinv*qinv*qinv*qinv  + fourmassSquered * (qinv*qinv - r*r );
 571   if (temp < 0.0)
 572    {
 573      Error("GetOneD","temp is less then 0: %f",temp);
 574      return;
 575    }
 576   temp = temp*(msqrd+firstPx*firstPx);
 577
 578   px = (px - TMath::Sqrt(temp))/(2.*msqrd);
 579
 580   Double_t py = r*TMath::Sin(phi);
 581   Double_t pz = r*TMath::Cos(phi);
 582
 583   TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
 584   TVector3 vector(px,py,pz);
 585   Rotate(firstpvector,vector);
 586
 587   Double_t e = TMath::Sqrt(msqrd + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
 588   second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
 589 //  TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, firstPx,0,0,e=TMath::Sqrt(msqrd+firstPx*firstPx),0.0,0.0,0.0,0.0);
 590 //        TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
 591 //                px, py, pz, e, vx, vy, vz, tof);
 592
 593   if (GetDebug()) Info("GetOneD","Randomized qinv = %f, obtained = %f",qinv,GetQInv(first,second));
 594
 595 }
 596 /***********************************************************/
 597
 598 Int_t AliGenHBTosl::GetThreeD(TParticle* first,TParticle* second, Double_t qout, Double_t qside, Double_t qlong)
 599 {
 600 //deprecated method that caclulates momenta of the second particle
 601 //out of  qout qside and qlong and the first particle
 602   Double_t m = first->GetMass();
 603   Double_t m2 = m*m;
 604
 605   Double_t px = first->P();//first particle is rotated that only X is non-zero  thus P==Px
 606   Double_t px2 = px*px;
 607
 608
 609   Double_t qout2 = qout*qout;
 610   Double_t qside2 = qside*qside;
 611   Double_t qlong2 = qlong*qlong;
 612
 613
 614   Double_t util1 = 4.*px2 - qside2;//4*P1x^2 - Y^2
 615   if (util1 < 0)
 616    {
 617      Info("GetThreeD","4.*px2* - qside2 is negative px: %f, qside: %f",px,qside);
 618      return 1;
 619    }
 620   Double_t util2 = TMath::Sqrt(px2*qout2*util1);
 621
 622
 623   Double_t p2x,p2y,p2z;
 624
 625 //  if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))
 626   if (qout >= 0)
 627    {
 628      //p2x
 629      Double_t tmp = px*(2.*px2 - qside2);
 630      tmp -=  util2;
 631      p2x = tmp/(2.*px2);
 632
 633      //p2y
 634      tmp =  qout - TMath::Sqrt(util1);
 635      p2y = - (tmp*qside)/(2.*px);
 636
 637      //p2z
 638      tmp = 4.*m2 + 2.*qout2+qlong2;
 639      tmp *= px;
 640      tmp -= 2.*util2;//!!!
 641      tmp += 4.*px*px2;
 642      tmp *= qlong2;
 643
 644      Double_t m2px2 = m2+px2;
 645      Double_t tmp2 = m2px2*tmp;
 646
 647      tmp  = 4.*(m2px2+qout2) + qlong2;
 648      tmp *= px;
 649      tmp -= 4.*util2;
 650      tmp *= 4.*(m2px2) + qlong2;
 651      tmp *= qlong2*qlong2;
 652      tmp *= m2px2*m2px2;
 653      tmp *= px;
 654      if (tmp < 0)
 655       {
 656         Error("","Argument of sqrt is negative");
 657         return 1;
 658       }
 659
 660      tmp2 += TMath::Sqrt(tmp);
 661
 662      tmp = 8.0*px*m2px2*m2px2;
 663      p2z = -TMath::Sqrt(tmp2/tmp);
 664      if (qlong < 0) p2z = -p2z;
 665    }
 666   else
 667    {
 668      //p2x
 669      Double_t tmp = px*(2.*px2 - qside2);
 670      tmp +=  util2;
 671      p2x = tmp/(2.*px2);
 672
 673      //p2y
 674      tmp =  qout - TMath::Sqrt(util1);
 675      p2y = - (tmp*qside)/(2.*px);
 676
 677      //p2z
 678      tmp = 4.*m2 + 2.*qout2+qlong2;
 679      tmp *= px;
 680      tmp += 2.*util2;//!!!
 681      tmp += 4.*px*px2;
 682      tmp *= qlong2;
 683
 684      Double_t m2px2 = m2+px2;
 685      Double_t tmp2 = m2px2*tmp;
 686
 687      tmp  = 4.*(m2px2+qout2) + qlong2;
 688      tmp *= px;
 689      tmp += 4.*util2;
 690      tmp *= 4.*(m2px2) + qlong2;
 691      tmp *= qlong2*qlong2;
 692      tmp *= m2px2*m2px2;
 693      tmp *= px;
 694      if (tmp < 0)
 695       {
 696         Error("","Argument of sqrt is negative");
 697         return 1;
 698       }
 699
 700      tmp2 += TMath::Sqrt(tmp);
 701
 702      tmp = 8.0*px*m2px2*m2px2;
 703      p2z = -TMath::Sqrt(tmp2/tmp);
 704      if (qlong < 0) p2z = -p2z;
 705     }
 706
 707 //  if ( (qside >= 0) && (qout >= 0) && (qlong >= 0))  p2z = -p2z;
 708
 709   TVector3 firstpvector(first->Px(),first->Py(),first->Pz());
 710   TVector3 vector(p2x,p2y,p2z);
 711   Rotate(firstpvector,vector);
 712
 713   Double_t e = TMath::Sqrt(m2 + vector.X()*vector.X() + vector.Y()*vector.Y() + vector.Z()*vector.Z());
 714   second->SetMomentum(vector.X(),vector.Y(),vector.Z(),e);
 715
 716 ////////////
 717   if ( GetDebug() > 3 )
 718    {
 719      e=TMath::Sqrt(m2+px*px);
 720      TParticle* f = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, px , 0.0, 0.0, e,0.0,0.0,0.0,0.0);
 721
 722      e = TMath::Sqrt(m2 + p2x*p2x + p2y*p2y + p2z*p2z);
 723      TParticle* s = new TParticle(first->GetPdgCode(),0,-1,-1,-1,-1, p2x, p2y, p2z, e, 0.0, 0.0, 0.0, 0.0);
 724
 725      Double_t qo, qs, ql;
 726      GetQOutQSideQLong(f,s,qo, qs, ql);
 727
 728      Info("GetThreeD","TEST");
 729      f->Print();
 730      s->Print();
 731      Info("GetThreeD","Required %f %f %f",qout,qside,qlong);
 732      Info("GetThreeD","Got      %f %f %f",qo,qs,ql);
 733      if ( qout == qo)
 734         if (qside == qs)
 735           if  (qlong == ql)
 736             Info("GetThreeD","!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
 737    }
 738 /////////////
 739   return 0;
 740 }
 741 /***********************************************************/
 742
 743 void AliGenHBTosl::StartSignal()
 744 {
 745 //Starts the signal histograms
 746  ofstream& logfile = *fLogFile;
 747  logfile<<"\n\n\n\n";
 748  logfile<<"************************************************"<<endl;
 749  logfile<<"************************************************"<<endl;
 750  logfile<<"               StartSignal                      "<<endl;
 751  logfile<<"************************************************"<<endl;
 752  logfile<<"************************************************"<<endl;
 753
 754  AliStack* stack;
 755
 756  fSwapped = kFALSE;
 757
 758  TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
 759  TParticle* second = &particle;
 760
 761  TIter next(fStackBuffer);
 762  while(( stack=(AliStack*)next() ))
 763   {
 764     stack->Reset();
 765   }
 766
 767  AliStack* genstack = fGenerator->GetStack();
 768  if (genstack == 0x0)
 769   {
 770    genstack = new AliStack(fNpart);
 771    fGenerator->SetStack(genstack);
 772   }
 773  else
 774   {
 775    genstack->Reset();
 776   }
 777
 778  StartSignalPass1();
 779  //We alread have detailed histograms and we do not need Coarse anymore
 780  delete fQCoarseSignal;
 781  delete fQCoarseBackground;
 782  fQCoarseSignal = 0x0;
 783  fQCoarseBackground = 0x0;
 784
 785
 786  const Double_t kNDF = fQNBins*fQNBins*fQNBins;
 787
 788  TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
 789  work->Sumw2();
 790  work->SetDirectory(0x0);
 791
 792
 793  Double_t binwdh = work->GetBinWidth(1)/2.;
 794
 795  Double_t*** chiarray = new Double_t** [fQNBins+1];
 796  Double_t*** sigarray = new Double_t** [fQNBins+1];
 797
 798  for (Int_t i = 1; i<=fQNBins; i++)
 799    {
 800       chiarray[i] =  new Double_t* [fQNBins+1];
 801       sigarray[i] =  new Double_t* [fQNBins+1];
 802
 803       for (Int_t k = 1; k<=fQNBins; k++)
 804        {
 805          chiarray[i][k] =  new Double_t [fQNBins+1];
 806          sigarray[i][k] =  new Double_t [fQNBins+1];
 807        }
 808    }
 809
 810
 811  Float_t chisqrchange = fMaxChiSquereChange + 1.;
 812  Float_t chisqrPerDF = fMaxChiSquerePerNDF;
 813  Float_t chisqrold = 0.0;
 814
 815  Int_t counter = 1;
 816  Int_t niterations = 1;
 817  Int_t rotaxisorder = 1;//defines order of looping over 3D histo (X,Y,Z or Y,Z,X or Z,X,Y)
 818
 819  Bool_t flag = kTRUE;
 820  Bool_t shortloop = kTRUE;
 821  TCanvas* c1 = new TCanvas();
 822
 823
 824  printf("\n");
 825  Info("StartSignal","\n\n\n\nSecond Pass\n\n\n\n");
 826
 827  while ( ( (chisqrPerDF > fMaxChiSquereChange) || flag) && (niterations++ < fMaxIterations)  )
 828   {
 829
 830    logfile<<"StartSignal\n";
 831    logfile<<" Row 1 Theory, 2 current value, 3 Chi2 \n";
 832
 833    Double_t chisqrnew = 0.0;
 834    flag = kFALSE;
 835
 836    Double_t scale = Scale(fQSignal,fQBackground);
 837    work->Divide(fQSignal,fQBackground,scale);
 838
 839    if ( (counter%100) == 0)
 840     {
 841       c1->cd();
 842       char buff[50];
 843       sprintf(buff,"QTWorkPass2.%3d.root",counter);
 844       TFile* file = TFile::Open(buff,"update");
 845       work->Write();
 846       work->SetDirectory(0x0);
 847       file->Close();
 848       delete file;
 849
 850       sprintf(buff,"QTBackgroundPass2.%3d.root",counter);
 851       file = TFile::Open(buff,"update");
 852       fQBackground->Write();
 853       fQBackground->SetDirectory(0x0);
 854       file->Close();
 855       delete file;
 856
 857       sprintf(buff,"QTSignalPass2.%3d.root",counter);
 858       file = TFile::Open(buff,"update");
 859       fQSignal->Write();
 860       fQSignal->SetDirectory(0x0);
 861       file->Close();
 862       delete file;
 863     }
 864
 865    counter++;
 866    Int_t novertresh = 0;
 867    for (Int_t k = 1; k<=fQNBins; k++)
 868     {
 869       Double_t z = work->GetZaxis()->GetBinCenter(k);
 870       for (Int_t j = 1; j<=fQNBins; j++)
 871         {
 872           Double_t y = work->GetYaxis()->GetBinCenter(j);
 873           for (Int_t i = 1; i<=fQNBins; i++)
 874             {
 875               Double_t x = work->GetXaxis()->GetBinCenter(i);//get center value of a bin (qout)
 876               sigarray[i][j][k] = fQSignal->GetBinContent(i,j,k);//store current value of signal histogram
 877               Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);//get expected value of CF in that qinv
 878               Double_t diff = v - work->GetBinContent(i,j,k);//store difference betweeon current value, and desired value
 879               chiarray[i][j][k] = diff; // no-x x is a weight to get good distribution
 880               Double_t be = work->GetBinError(i,j,k);
 881               chisqrnew += diff*diff/(be*be);//add up chisq
 882
 883               //even if algorithm is stable (chi sqr change less then threshold)
 884               //and any bin differs more then 5% from expected value we continue
 885               Double_t fact = diff;
 886               if (TMath::Abs(fact) > 0.1)
 887                {
 888                  flag = kTRUE;
 889                  novertresh++;
 890                }
 891             }
 892          }
 893      }
 894
 895
 896    char msg[1000];
 897
 898    printf("\n");
 899
 900    for (Int_t k = 25; k < 36; k++)
 901     {
 902       Double_t tx = work->GetXaxis()->GetBinCenter(30);
 903       Double_t ty = work->GetYaxis()->GetBinCenter(30);
 904       Double_t tz = work->GetZaxis()->GetBinCenter(k);
 905       sprintf(msg,"% 6.5f ",GetQOutQSideQLongCorrTheorValue(tx,ty,tz));
 906       logfile<<msg;
 907     }
 908    logfile<<endl;
 909
 910    for (Int_t k = 25; k < 36; k++)
 911     {
 912       sprintf(msg,"%6.5f ",work->GetBinContent(30,30,k));
 913       logfile<<msg;
 914     }
 915    logfile<<endl;
 916
 917    for (Int_t k = 25; k < 36; k++)
 918     {
 919       sprintf(msg,"% 6.5f ",chiarray[30][30][k]);
 920       logfile<<msg;
 921     }
 922    logfile<<endl;
 923
 924    chisqrchange = TMath::Abs(chisqrnew - chisqrold)/chisqrnew;
 925    chisqrold = chisqrnew;
 926
 927    chisqrPerDF = chisqrnew/kNDF;
 928
 929    Info("StartSignal","Iteration %d Chi-sq change %f%%",niterations,chisqrchange*100.0);
 930    Info("StartSignal","ChiSq = %f, NDF = %f, ChiSq/NDF = %f",chisqrnew, kNDF, chisqrPerDF );
 931    Info("StartSignal","novertresh = %d",novertresh);
 932
 933
 934    stack = RotateStack();
 935    genstack->Reset();
 936    fGenerator->Generate();
 937    Int_t ninputparticle = 0, ntr = 0;
 938    if ( genstack->GetNtrack() < fNpart/2)
 939     {
 940       Warning("StartSignal","**********************************");
 941       Warning("StartSignal","Generator generated (%d) less particles then expected (%d).",
 942                genstack->GetNtrack(),fNpart/2);
 943       Warning("StartSignal","**********************************");
 944     }
 945
 946    Int_t sc = 0; //security check against infinite loop
 947    while ( (ntr+1) < fNpart)//ntr is number of track generated up to now
 948     {
 949      Int_t xmax = 1;
 950      Int_t ymax = 1;
 951      Int_t zmax = 1;
 952      Double_t qout;
 953      Double_t qside;
 954      Double_t qlong;
 955
 956      Int_t i,j,k;
 957
 958      Int_t* cx = 0x0;
 959      Int_t* cy = 0x0;
 960      Int_t* cz = 0x0;
 961
 962      switch (rotaxisorder)
 963       {
 964         case 1:
 965           cx = &i;
 966           cy = &j;
 967           cz = &k;
 968           break;
 969         case 2:
 970           cx = &j;
 971           cy = &k;
 972           cz = &i;
 973           break;
 974         case 3:
 975           cx = &k;
 976           cy = &i;
 977           cz = &j;
 978           break;
 979       }
 980      rotaxisorder++;
 981      if (rotaxisorder > 3) rotaxisorder = 1;
 982      Int_t nrange;
 983
 984      if (shortloop)
 985       {
 986         shortloop = kFALSE;
 987         nrange = fQNBins;
 988       }
 989      else
 990       {
 991         shortloop = kTRUE;
 992         nrange = fQNBins/4;
 993       }
 994
 995 //     Bool_t force = kFALSE;
 996      for ( k = 1; k <=nrange;k++ )
 997       {
 998        for ( j = 1; j<=nrange; j++)
 999          {
1000           for ( i = 1; i<=nrange; i++)
1001             {
1002               if ( (chiarray[*cx][*cy][*cz]) > (chiarray[xmax][ymax][zmax]) )
1003                {
1004                  xmax = *cx;
1005                  ymax = *cy;
1006                  zmax = *cz;
1007                }
1008
1009 //              Double_t fact = chiarray[*cx][*cy][*cz];//chiarray is chi2*qinv
1010 //              if (fact > work->GetBinError(*cx,*cy,*cz))//if differece between what we want and
1011 //               {                                  //what we have is bigger than stat. error
1012 //                                                  //we force to fill that bin
1013 //                 qout  = work->GetXaxis()->GetBinCenter(*cx);
1014 //                 qside = work->GetYaxis()->GetBinCenter(*cy);
1015 //                 qlong = work->GetZaxis()->GetBinCenter(*cz);
1016
1017 //                 Info("StartSignal"," bin: (%d,%d,%d) loop status (%d,%d,%d) \nUsing Force: chiarray: %f \nq(o,s,l): (%f,%f,%f)  signal: %d background: %d binerror: %f",
1018 //        *cx,*cy,*cz,i,j,k,fact,qout,qside,qlong,
1019 //        (Int_t)sigarray[*cx][*cy][*cz],(Int_t)fQBackground->GetBinContent(*cx,*cy,*cz),work->GetBinError(*cx,*cy,*cz));
1020 //                 force = kTRUE;
1021 //                 break;
1022 //               }
1023
1024             }
1025 //           if (force) break;
1026           }
1027 //         if (force) break;
1028         }
1029
1030
1031       qout  = work->GetXaxis()->GetBinCenter(xmax);
1032       qside = work->GetYaxis()->GetBinCenter(ymax);
1033       qlong = work->GetZaxis()->GetBinCenter(zmax);
1034
1035 //      Info("StartSignal"," bin: (%d,%d,%d) chiarray: %f \nq(o,s,l): (%f,%f,%f)  signal: %d background: %d binerror: %f",
1036 //            xmax,ymax,zmax,chiarray[xmax][ymax][zmax],qout,qside,qlong,
1037 //            (Int_t)sigarray[xmax][ymax][zmax],
1038 //            (Int_t)fQBackground->GetBinContent(xmax,ymax,zmax),
1039 //            work->GetBinError(xmax,ymax,zmax));
1040
1041       qout  = gRandom->Uniform(qout-binwdh,qout+binwdh);
1042       qside = gRandom->Uniform(qside-binwdh,qside+binwdh);
1043       qlong = gRandom->Uniform(qlong-binwdh,qlong+binwdh);
1044
1045       TParticle* first = 0;
1046       while (ninputparticle < genstack->GetNtrack())
1047        {
1048          TParticle* tmpp = genstack->Particle(ninputparticle++);
1049          if (tmpp->GetPdgCode() == fPID)
1050           {
1051            if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1052             {
1053               first = tmpp;
1054               break;
1055             }
1056           }
1057        }
1058
1059       if (first == 0x0)
1060        {
1061          if ( fDebug > 2 ) Info("StartSignal","No more particles of that type");
1062          break;
1063        }
1064
1065       Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1066       if (retval)
1067        {
1068          Info("StartSignal","Can not find momenta for this OSL and particle OSL = %f %f %f",qout,qside,qlong);
1069          first->Print();
1070          second->Print();
1071
1072          continue;
1073        }
1074      //in case this particle is falling into signal area with another
1075      //particle we take a another pair
1076      //it can intruduce artificial correlations
1077       Bool_t checkresult = CheckParticle(second,first,stack);
1078       if ( checkresult  && (sc < 10) )
1079        {
1080          sc++;
1081          continue;
1082        }
1083       sc = 0;
1084
1085       //Put on output stack
1086       SetTrack(first,ntr,stack);
1087       SetTrack(second,ntr,stack);
1088
1089       Double_t y = GetQOutQSideQLongCorrTheorValue(qout,qside,qlong);
1090
1091       sigarray[xmax][ymax][zmax] ++;
1092       chiarray[xmax][ymax][zmax] = scale*sigarray[xmax][ymax][zmax]/fQBackground->GetBinContent(xmax,ymax,zmax);
1093       chiarray[xmax][ymax][zmax] = (y - chiarray[xmax][ymax][zmax]);
1094
1095     }
1096    Info("StartSignal","Mixing background...");
1097    Mix(fStackBuffer,fQBackground,fQSecondBackground); //upgrate background
1098    Info("StartSignal","Mixing signal...");
1099    Mix(stack,fQSignal,fQSecondSignal); //upgrate background
1100
1101
1102 //   if ( (chisqrPerDF < 2.0) && (fSwapped == kFALSE) )
1103 //     {
1104 //         SwapGeneratingHistograms();
1105 //     }
1106
1107  }
1108  TFile* filef = TFile::Open("QTBackground.root","recreate");
1109  fQBackground->Write();
1110  fQBackground->SetDirectory(0x0);
1111  filef->Close();
1112  delete filef;
1113
1114  filef = TFile::Open("QTSignal.root","recreate");
1115  fQSignal->Write();
1116  fQSignal->SetDirectory(0x0);
1117  filef->Close();
1118  delete filef;
1119
1120
1121  delete c1;
1122  delete work;
1123
1124  for (Int_t i = 1; i<=fQNBins; i++)
1125    {
1126      for (Int_t k = 1; k<=fQNBins; k++)
1127       {
1128         delete [] chiarray[i][k];
1129         delete [] sigarray[i][k];
1130       }
1131      delete [] chiarray[i];
1132      delete [] sigarray[i];
1133    }
1134  delete [] chiarray;
1135  delete [] sigarray;
1136
1137 }
1138 /***********************************************************/
1139
1140 void AliGenHBTosl::StartSignalPass1()
1141 {
1142  //This method makes first part of the initialization of working histograms
1143  //It randomizes qout, qside and qlong from the coarse signal histogram
1144
1145  Bool_t flag = kTRUE;
1146  TParticle particle(fPID,0,-1,-1,-1,-1,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
1147  TParticle* second = &particle;
1148  Double_t qout;
1149  Double_t qside;
1150  Double_t qlong;
1151
1152  Info("StartSignalPass1","\n\nFirst Pass\n\n");
1153
1154  while (flag)
1155   {
1156     Info("StartSignalPass1","NextEvent");
1157     AliStack* stack = RotateStack();
1158     AliStack* genstack = fGenerator->GetStack();
1159     genstack->Reset();
1160     fGenerator->Generate();
1161     Int_t j = 0, ntr = 0;
1162     if ( genstack->GetNtrack() < fNpart/2)
1163      {
1164        Warning("StartSignalPass1","**********************************");
1165        Warning("StartSignalPass1","Generator generated (%d) less particles then expected (%d).",
1166                 genstack->GetNtrack(),fNpart/2);
1167        Warning("StartSignalPass1","**********************************");
1168      }
1169
1170     Int_t sc = 0;//security check against infinite loop
1171     while ((ntr+1)<fNpart)
1172      {
1173
1174 //    Info("StartSignal","Number of track on output stack: = %d", ntr);
1175 //    Info("StartSignal","Number of track on input stack: = %d\n", j);
1176
1177       TParticle* first = 0;
1178       while (j < genstack->GetNtrack())
1179        {
1180          TParticle* tmpp = genstack->Particle(j++);
1181          if (tmpp->GetPdgCode() == fPID)
1182           {
1183            if (CheckParticle(tmpp,0x0,stack) == kFALSE)
1184             {
1185               first = tmpp;
1186               break;
1187             }
1188            else
1189             {
1190               Info("StartSignalPass1","Particle did not pass the safety check 1");
1191               tmpp->Print();
1192             }
1193           }
1194        }
1195
1196       if (first == 0x0)
1197        {
1198          if ( fDebug > 2 ) Info("StartSignalPass1","No more particles of that type");
1199
1200          break;
1201        }
1202
1203       SetTrack(first,ntr,stack);
1204
1205       fQCoarseSignal->GetRandom3(qout,qside,qlong);
1206
1207       Int_t retval = GetThreeD(first,second,qout,qside,qlong);
1208       if (retval)
1209        {
1210          //Info("StartSignal","Can not find momenta for this OSL and particle");
1211          continue;
1212        }
1213       //in case this particle is falling into signal area with another
1214       //particle we take a another pair
1215       //it can intruduce artificial correlations
1216        Bool_t checkresult = CheckParticle(second,first,stack);
1217        if ( checkresult  && (sc < 10) )
1218         {
1219           sc++;
1220           Info("StartSignalPass1","Particle did not pass the safety check 2");
1221           second->Print();
1222           continue;
1223         }
1224
1225        sc = 0;
1226
1227       SetTrack(second,ntr,stack);
1228      }
1229
1230     Mix(stack,fQSignal,fQSecondSignal);
1231     Mix(fStackBuffer,fQBackground,fQSecondBackground);
1232
1233     flag = kFALSE;
1234
1235     for (Int_t k = 1; k<=fQNBins; k++)
1236       {
1237        for (Int_t j = 1; j<=fQNBins; j++)
1238          {
1239            for (Int_t i = 1; i<=fQNBins; i++)
1240              {
1241                 if ( (fQBackground->GetBinContent(i,j,k) < fMinFill) )
1242                   {
1243                     //(fQSignal->GetBinContent(i,j,k) < fMinFill) ||
1244         Info("StartSignalPass1","bin (%d,%d,%d): signal=%f background=%f",i,j,k,
1245                           fQSignal->GetBinContent(i,j,k),fQBackground->GetBinContent(i,j,k));
1246         flag = kTRUE;//continue while
1247         break;//breakes for not while
1248                   }
1249              }
1250             if (flag == kTRUE) break;
1251           }
1252         if (flag == kTRUE) break;
1253       }
1254
1255   }//while (flag)
1256
1257
1258 }
1259 /***********************************************************/
1260
1261 void AliGenHBTosl::FillCoarseSignal()
1262 {
1263  //Makes coarse signal by multiplying the coarse background and required function
1264  Info("FillCoarseSignal","START");
1265  for (Int_t k = 1; k <=fQNBins ;k++ )
1266    {
1267     Double_t z = fQCoarseBackground->GetZaxis()->GetBinCenter(k);
1268     for (Int_t j = 1; j <=fQNBins; j++)
1269       {
1270        Double_t y = fQCoarseBackground->GetYaxis()->GetBinCenter(j);
1271        for (Int_t i = 1; i <=fQNBins; i++)
1272          {
1273            Double_t x  = fQCoarseBackground->GetXaxis()->GetBinCenter(i);
1274            Double_t v = GetQOutQSideQLongCorrTheorValue(x,y,z);
1275            Info("FillCoarseSignal","Bin (%d,%d,%d): osl(%f,%f,%f)=%f",i,j,k,x,y,z,v);
1276            fQCoarseSignal->SetBinContent(i,j,k,v*fQCoarseBackground->GetBinContent(i,j,k));
1277          }
1278        }
1279     }
1280
1281   //if (GetDebug())
1282   TestCoarseSignal();
1283
1284   Info("FillCoarseSignal","DONE");
1285 }
1286 /***********************************************************/
1287
1288 void AliGenHBTosl::FillCoarse()
1289 {
1290   //creates the statistical background histogram on the base of input from
1291   //fGenerator
1292   Info("FillCoarse","START");
1293
1294   AliStack* stack;
1295   Int_t niter = 0;
1296
1297   Bool_t cont;
1298   TH3D tmph("tmph","tmph",2,0,1,2,0,1,2,0,1);
1299   printf("\n");
1300
1301   do
1302    {
1303 //     if (niter > 20) break;
1304
1305      cout<<niter++<<"  bincont "<<fQCoarseBackground->GetBinContent(30,30,28)
1306                   <<" "<<fQCoarseBackground->GetBinContent(30,30,29)
1307                   <<" "<<fQCoarseBackground->GetBinContent(30,30,30)
1308                   <<" "<<fQCoarseBackground->GetBinContent(30,30,31)
1309                   <<" "<<fQCoarseBackground->GetBinContent(30,30,32)
1310                   <<"\n";
1311      fflush(0);
1312
1313      stack = RotateStack();
1314      fGenerator->SetStack(stack);
1315      fGenerator->Init();
1316      fGenerator->Generate();
1317
1318      Mix(fStackBuffer,fQCoarseBackground,&tmph);
1319
1320      cont = kFALSE;
1321
1322      Info("FillCoarse","fMinFill = %d",fMinFill);
1323      for (Int_t k = 1; k<=fQNBins; k++)
1324        {
1325         for (Int_t j = 1; j<=fQNBins; j++)
1326           {
1327             for (Int_t i = 1; i<=fQNBins; i++)
1328               {
1329                 if ( fQCoarseBackground->GetBinContent(i,j,k) < fMinFill )
1330                 {
1331                   cont = kTRUE;
1332                   Info("FillCoarse","bin (%d,%d,%d)=%f",i,j,k,fQCoarseBackground->GetBinContent(i,j,k));
1333                   break;
1334                 }
1335
1336               }
1337             if (cont) break;
1338           }
1339          if (cont) break;
1340         }
1341    }while(cont);
1342
1343   printf("\n");
1344   fGenerator->SetStack(0x0);
1345   Info("FillCoarse","DONE");
1346
1347 }
1348 /***********************************************************/
1349
1350 void AliGenHBTosl::Mix(TList* eventbuffer,TH3D* denominator,TH3D* denominator2)
1351 {
1352   //Fills denominators
1353   //Mixes events stored in the eventbuffer and fills the background histograms
1354   static TStopwatch stoper;
1355
1356   if (eventbuffer == 0x0)
1357    {
1358     Error("Mix","Buffer List is null.");
1359     return;
1360    }
1361
1362   if (denominator == 0x0)
1363    {
1364     Error("Mix","Denominator histogram is null.");
1365     return;
1366    }
1367
1368   if (denominator2 == 0x0)
1369    {
1370     Error("Mix","Denominator2 histogram is null.");
1371     return;
1372    }
1373
1374   Info("Mix","%s",denominator->GetName());
1375   stoper.Start();
1376
1377   TIter next(eventbuffer);
1378   AliStack* firstevent;
1379   AliStack* secondevent = 0x0;
1380
1381   while(( firstevent=(AliStack*)next() ))
1382    {
1383     if (secondevent == 0x0)
1384      {
1385        secondevent = firstevent;
1386        continue;
1387      }
1388 //    Info("Mix","Mixing %#x with %#x",firstevent,secondevent);
1389     for(Int_t j = 0; j < firstevent->GetNtrack(); j++ )
1390      {
1391        TParticle* firstpart = firstevent->Particle(j);
1392
1393        Float_t phi = firstpart->Phi();
1394        if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1395
1396 //       Info("Mix","Mixing %d phi %f min %f max %f",j,phi,fSamplePhiMin,fSamplePhiMax);
1397
1398        for(Int_t i = 0; i < secondevent->GetNtrack(); i++ )
1399          {
1400            TParticle* secondpart = secondevent->Particle(i);
1401            phi = secondpart->Phi();
1402            if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1403
1404            Double_t qout;
1405            Double_t qside;
1406            Double_t qlong;
1407            GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1408            denominator->Fill(qout,qside,qlong);
1409            denominator2->Fill(qout,qside,qlong);
1410          }
1411      }
1412
1413     secondevent = firstevent;
1414    }
1415   stoper.Stop();
1416   stoper.Print();
1417
1418 }
1419 /***********************************************************/
1420
1421 void AliGenHBTosl::Mix(AliStack* stack, TH3D* numerator, TH3D* numerator2)
1422 {
1423 //fils numerator with particles from stack
1424   static TStopwatch stoper;
1425   if (stack == 0x0)
1426    {
1427     Error("Mix","Stack is null.");
1428     return;
1429    }
1430
1431   if ( (numerator == 0x0) || (numerator2 == 0x0) )
1432    {
1433     Error("Mix","Numerator histogram is null.");
1434     return;
1435    }
1436
1437   Info("Mix","%s",numerator->GetName());
1438   stoper.Start();
1439
1440   for(Int_t j = 0; j < stack->GetNtrack(); j++ )
1441     {
1442       TParticle* firstpart = stack->Particle(j);
1443       Float_t phi = firstpart->Phi();
1444       if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1445
1446       for(Int_t i = j+1; i < stack->GetNtrack(); i++ )
1447        {
1448          TParticle* secondpart = stack->Particle(i);
1449          phi = secondpart->Phi();
1450          if ( (phi < fSamplePhiMin) || ( phi > fSamplePhiMax)) continue;
1451          Double_t qout;
1452          Double_t qside;
1453          Double_t qlong;
1454          GetQOutQSideQLong(firstpart,secondpart,qout,qside,qlong);
1455          numerator->Fill(qout,qside,qlong);
1456          numerator2->Fill(qout,qside,qlong);
1457        }
1458     }
1459   stoper.Stop();
1460   stoper.Print();
1461
1462 }
1463 /***********************************************************/
1464
1465 Double_t AliGenHBTosl::GetQInv(TParticle* f, TParticle* s)
1466 {
1467 //calculates qinv
1468 // cout<<f->Px()<<"   "<<s->Px()<<endl;
1469  Double_t pxdiff = f->Px() - s->Px();
1470  Double_t pydiff = f->Py() - s->Py();
1471  Double_t pzdiff = f->Pz() - s->Pz();
1472  Double_t ediff  = f->Energy() - s->Energy();
1473
1474  Double_t qinvl = ediff*ediff - ( pxdiff*pxdiff + pydiff*pydiff + pzdiff*pzdiff );
1475  Double_t qinv = TMath::Sqrt(TMath::Abs(qinvl));
1476  return qinv;
1477 }
1478 /***********************************************************/
1479
1480 void  AliGenHBTosl::GetQOutQSideQLong(TParticle* f, TParticle* s,Double_t& out, Double_t& side, Double_t& lon)
1481 {
1482  //returns qout,qside and qlong of the pair of particles
1483  out = side = lon = 10e5;
1484
1485  Double_t pxsum = f->Px() + s->Px();
1486  Double_t pysum = f->Py() + s->Py();
1487  Double_t pzsum = f->Pz() + s->Pz();
1488  Double_t esum  = f->Energy() + s->Energy();
1489  Double_t pxdiff = f->Px() - s->Px();
1490  Double_t pydiff = f->Py() - s->Py();
1491  Double_t pzdiff = f->Pz() - s->Pz();
1492  Double_t ediff  = f->Energy() - s->Energy();
1493  Double_t kt =  0.5*TMath::Hypot(pxsum,pysum);
1494
1495  Double_t k2 = pxsum*pxdiff+pysum*pydiff;
1496
1497  if (kt == 0.0)
1498   {
1499      f->Print();
1500      s->Print();
1501      kt = 10e5;
1502   }
1503  else
1504   {
1505     out = 0.5*k2/kt;
1506     side = (f->Px()*s->Py()-s->Px()*f->Py())/kt;
1507   }
1508
1509  Double_t beta = pzsum/esum;
1510  Double_t gamma = 1.0/TMath::Sqrt(1.0 - beta*beta);
1511
1512  lon = gamma * ( pzdiff - beta*ediff );
1513
1514 // out = TMath::Abs(out);
1515 // side = TMath::Abs(side);
1516 // lon = TMath::Abs(lon);
1517 }
1518
1519 /***********************************************************/
1520
1521 Double_t AliGenHBTosl::Scale(TH3D* num, TH3D* den)
1522 {
1523  //Calculates the factor that should be used to scale
1524  //quatience of num and den to 1 at tail
1525
1526   if (GetDebug()) Info("Scale","Entered Scale()");
1527   if(!num)
1528    {
1529      Error("Scale","No numerator");
1530      return 0.0;
1531    }
1532   if(!den)
1533    {
1534      Error("Scale","No denominator");
1535      return 0.0;
1536    }
1537
1538   if(fNBinsToScale < 1)
1539    {
1540     return 0.0;
1541     Error("Scale","Number of bins for scaling is smaller thnan 1");
1542    }
1543   Int_t fNBinsToScaleX = fNBinsToScale;
1544   Int_t fNBinsToScaleY = fNBinsToScale;
1545   Int_t fNBinsToScaleZ = fNBinsToScale;
1546
1547   Int_t nbinsX = num->GetNbinsX();
1548   if (fNBinsToScaleX > nbinsX)
1549    {
1550     Error("Scale","Number of X bins for scaling is bigger thnan number of bins in histograms");
1551     return 0.0;
1552    }
1553
1554   Int_t nbinsY = num->GetNbinsX();
1555   if (fNBinsToScaleY > nbinsY)
1556    {
1557     Error("Scale","Number of Y bins for scaling is bigger thnan number of bins in histograms");
1558     return 0.0;
1559    }
1560
1561   Int_t nbinsZ = num->GetNbinsZ();
1562   if (fNBinsToScaleZ > nbinsZ)
1563    {
1564     Error("Scale","Number of Z bins for scaling is bigger thnan number of bins in histograms");
1565     return 0.0;
1566    }
1567
1568   if (GetDebug()>0) Info("Scale","No errors detected");
1569
1570   Int_t offsetX = nbinsX - fNBinsToScaleX - 1; //bin that we start loop over bins in axis X
1571   Int_t offsetY = nbinsY - fNBinsToScaleY - 1; //bin that we start loop over bins in axis Y
1572   Int_t offsetZ = nbinsZ - fNBinsToScaleZ - 1; //bin that we start loop over bins in axis Z
1573
1574   Double_t densum = 0.0;
1575   Double_t numsum = 0.0;
1576
1577   for (Int_t k = offsetZ; k<nbinsZ; k++)
1578     for (Int_t j = offsetY; j<nbinsY; j++)
1579       for (Int_t i = offsetX; i<nbinsX; i++)
1580        {
1581         if ( num->GetBinContent(i,j,k) > 0.0 )
1582          {
1583
1584            densum += den->GetBinContent(i,j,k);
1585            numsum += num->GetBinContent(i,j,k);
1586          }
1587        }
1588
1589   if(GetDebug() > 0)
1590     Info("Scale","numsum=%f densum=%f fNBinsToScaleX=%d fNBinsToScaleY=%d fNBinsToScaleZ=%d",
1591           numsum,densum,fNBinsToScaleX,fNBinsToScaleY,fNBinsToScaleZ);
1592
1593   if (numsum == 0) return 0.0;
1594   Double_t ret = densum/numsum;
1595
1596   if(GetDebug() > 0) Info("Scale","returning %f",ret);
1597   return ret;
1598
1599 }
1600 /***********************************************************/
1601
1602 void AliGenHBTosl::TestCoarseSignal()
1603 {
1604 //Tests how works filling from generated histogram shape
1605   TH3D* work = new TH3D("work","work",fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange,fQNBins,-fQRange,fQRange);
1606
1607 //  for (Int_t i = 0; i < fQCoarseBackground->GetEntries() ;i++)
1608 //   {
1609 //     Double_t x,y,z;
1610 //     fQCoarseSignal->GetRandom3(x,y,z);
1611 //     work->Fill(x,y,z);
1612 //   }
1613
1614   TCanvas* c1 = new TCanvas();
1615   c1->cd();
1616   work->Draw();
1617   c1->SaveAs("QTwork.root");
1618   TFile* file = TFile::Open("QTwork.root","update");
1619 //  work->Write();
1620   work->SetDirectory(0x0);
1621   file->Close();
1622
1623   fQCoarseSignal->Draw();
1624   c1->SaveAs("QTCoarseSignal.root");
1625   file = TFile::Open("QTCoarseSignal.root","update");
1626   fQCoarseSignal->Write();
1627   fQCoarseSignal->SetDirectory(0x0);
1628   file->Close();
1629
1630   fQCoarseBackground->Draw();
1631   c1->SaveAs("QTCoarseBackground.root");
1632   file = TFile::Open("QTCoarseBackground.root","update");
1633   fQCoarseBackground->Write();
1634   fQCoarseBackground->SetDirectory(0x0);
1635   file->Close();
1636
1637   TH1 *result = (TH1*)fQCoarseBackground->Clone("ratio");
1638   result->SetTitle("ratio");
1639   Float_t normfactor = Scale(work,fQCoarseBackground);
1640   result->Divide(work,fQCoarseBackground,normfactor);//work
1641
1642
1643   c1->cd();
1644   result->Draw();
1645   c1->SaveAs("QTresult.root");
1646   file = TFile::Open("QTresult.root","update");
1647   result->Write();
1648   result->SetDirectory(0x0);
1649   file->Close();
1650
1651   delete work;
1652   delete c1;
1653 }
1654 /***********************************************************/
1655
1656 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr)
1657 {
1658 //Shortcut to PushTrack(bla,bla,bla,bla.............)
1659    if (p->P() == 0.0)
1660     {
1661       Error("SetTrack(TParticle*,Int_t&)","Particle has zero momentum");
1662       return;
1663     }
1664
1665
1666    Int_t pdg = p->GetPdgCode();
1667    Double_t px = p->Px();
1668    Double_t py = p->Py();
1669    Double_t pz = p->Pz();
1670    Double_t e  = p->Energy();
1671    Double_t vx = p->Vx();
1672    Double_t vy = p->Vy();
1673    Double_t vz = p->Vz();
1674    Double_t tof = p->T();
1675
1676    TVector3 pol;
1677    p->GetPolarisation(pol);
1678
1679    Double_t polx = pol.X();
1680    Double_t poly = pol.Y();
1681    Double_t polz = pol.Z();
1682    TMCProcess mech = AliGenCocktailAfterBurner::IntToMCProcess(p->GetUniqueID());
1683    Float_t weight = p->GetWeight();
1684
1685    AliGenerator::PushTrack(fTrackIt, -1, pdg, px, py, pz, e, vx, vy, vz, tof,polx, poly, polz, mech, ntr, weight);
1686 }
1687 /***********************************************************/
1688
1689 void AliGenHBTosl::SetTrack(TParticle* p, Int_t& ntr, AliStack* stack) const
1690 {
1691 //Shortcut to SetTrack(bla,bla,bla,bla.............)
1692    if (p->P() == 0.0)
1693     {
1694       Error("SetTrack(TParticle*,Int_t&,AliStack*)","Particle has zero momentum");
1695       return;
1696     }
1697
1698    Int_t pdg = p->GetPdgCode();
1699    Double_t px = p->Px();
1700    Double_t py = p->Py();
1701    Double_t pz = p->Pz();
1702    Double_t e  = p->Energy();
1703
1704    stack->PushTrack(fTrackIt, -1, pdg, px, py, pz, e, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, kPPrimary, ntr,1,0);
1705 }
1706 /***********************************************************/
1707
1708 void AliGenHBTosl::Rotate(TVector3& relvector, TVector3& vector)
1709 {
1710 //This method rotates vector about the angeles that are needed to rotate
1711 //relvector from postion (firstPx,0,0) to its actual positon
1712 //In other words: To make equations easier
1713
1714   static TVector3 first;
1715   if (GetDebug())
1716    {
1717      first.SetXYZ(relvector.x(),relvector.y(),relvector.z());
1718    }
1719
1720   Double_t firstPx = TMath::Sqrt( relvector.x()*relvector.x() +
1721                                   relvector.y()*relvector.y() +
1722                                   relvector.z()*relvector.z()   );
1723
1724   Double_t rotAngleZ = -TMath::ATan2(relvector.y(),relvector.x());//calculating rot angles
1725   relvector.RotateZ(rotAngleZ);
1726   rotAngleZ = -rotAngleZ;
1727   Double_t rotAngleY = -TMath::ATan2(relvector.z(),relvector.x());
1728
1729   vector.RotateY(rotAngleY);
1730   vector.RotateZ(rotAngleZ);
1731
1732   if (GetDebug()>5)
1733    {
1734     TVector3 test(firstPx,0.0,0.0);
1735     test.RotateY(rotAngleY);
1736     test.RotateZ(rotAngleZ);
1737     ::Info("Rotate","Rotation test: px %f %f",first.x(),test.x());
1738     ::Info("Rotate","Rotation test: py %f %f",first.y(),test.y());
1739     ::Info("Rotate","Rotation test: pz %f %f",first.z(),test.z());
1740    }
1741 }
1742 /***********************************************************/
1743
1744 Double_t AliGenHBTosl::Rotate(Double_t x,Double_t y,Double_t z)
1745 {
1746 //Rotates vector to base where only x - coordinate is no-zero, and returns that
1747
1748   Double_t xylength = TMath::Hypot(x,y);
1749   Double_t sinphi = -y/xylength;
1750   Double_t cosphi = x/xylength;
1751
1752   Double_t xprime = cosphi*x - sinphi*y;
1753   Double_t yprime = sinphi*x + cosphi*y;
1754
1755   TVector3 v(x,y,z);
1756   Double_t a1 = -TMath::ATan2(v.Y(),v.X());
1757
1758   if (GetDebug()>5)
1759    {
1760      ::Info("Rotate","Xpr = %f  Ypr = %f",xprime,yprime);
1761      ::Info("Rotate","Calc sin = %f, and %f",sinphi,TMath::Sin(a1));
1762      ::Info("Rotate","Calc cos = %f, and %f",cosphi,TMath::Cos(a1));
1763    }
1764
1765   Double_t xprimezlength = TMath::Hypot(xprime,z);
1766
1767   Double_t sintheta = z/xprimezlength;
1768   Double_t costheta = xprime/xprimezlength;
1769
1770
1771   Double_t xbis = sintheta*z + costheta*(cosphi*x - sinphi*y);
1772
1773   ::Info("Rotate","Calculated rot %f, modulus %f",xbis,TMath::Sqrt(x*x+y*y+z*z));
1774   return xbis;
1775 }
1776 /***********************************************************/
1777
1778 AliStack* AliGenHBTosl::RotateStack()
1779 {
1780 //swaps to next stack last goes to first and is reseted
1781
1782  AliStack* stack;
1783  if ( fStackBuffer->GetSize() >= fBufferSize )
1784   {
1785     stack = (AliStack*)fStackBuffer->Remove(fStackBuffer->Last());
1786   }
1787  else
1788   {
1789     stack = new AliStack(fNpart);
1790   }
1791
1792  fStackBuffer->AddFirst(stack);
1793  stack->Reset();
1794  return stack;
1795 }
1796 /***********************************************************/
1797
1798 Double_t AliGenHBTosl::GetQInvCorrTheorValue(Double_t qinv) const
1799 {
1800 //Function (deprecated)
1801  static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1802
1803  return 1.0 + 0.5*TMath::Exp(-qinv*qinv*fQRadius*fQRadius/kFactorsqrd);
1804 }
1805 /***********************************************************/
1806
1807 Double_t AliGenHBTosl::GetQOutQSideQLongCorrTheorValue(Double_t& out, Double_t& side, Double_t& lon) const
1808 {
1809  //Theoretical function. Wa want to get correlation of the shape of this function
1810  static const Double_t kFactorsqrd = 0.197*0.197;//squared conversion factor SI<->eV
1811  return 1.0 + 0.7*TMath::Exp(-fQRadius*fQRadius*(out*out+side*side+lon*lon)/kFactorsqrd);
1812 }
1813 /***********************************************************/
1814
1815 Bool_t AliGenHBTosl::CheckParticle(TParticle* p, TParticle* aupair ,AliStack* stack)
1816 {
1817  //Checks if a given particle is falling into signal region with any other particle
1818  //already existing on stack
1819  return kFALSE;
1820
1821  if (fSignalRegion <=0) return kFALSE;
1822
1823  for (Int_t i = 0; i < stack->GetNtrack(); i++)
1824   {
1825     TParticle* part = stack->Particle(i);
1826     if (part == aupair) continue;
1827     Double_t qout = 10e5;
1828     Double_t qside= 10e5;
1829     Double_t qlong= 10e5;
1830     GetQOutQSideQLong(p,part,qout,qside,qlong);
1831
1832     if (TMath::Abs(qout)  < fSignalRegion)
1833       if (TMath::Abs(qside) < fSignalRegion)
1834        if (TMath::Abs(qlong) < fSignalRegion)
1835          return kTRUE;
1836   }
1837   return kFALSE;
1838 }
1839 /***********************************************************/
1840
1841 void AliGenHBTosl::SwapGeneratingHistograms()
1842 {
1843   //Checks if it is time to swap signal and background histograms
1844   //if yes it swaps them
1845   Int_t threshold = fMinFill;
1846   for (Int_t k = 1; k<=fQNBins; k++)
1847    {
1848      for (Int_t j = 1; j<=fQNBins; j++)
1849        {
1850          for (Int_t i = 1; i<=fQNBins; i++)
1851            {
1852              if ( fQSecondBackground->GetBinContent(i,j,k) < threshold) return;
1853            }
1854        }
1855
1856    }
1857
1858
1859   Info("SwapGeneratingHistograms","*******************************************");
1860   Info("SwapGeneratingHistograms","*******************************************");
1861   Info("SwapGeneratingHistograms","*******************************************");
1862   Info("SwapGeneratingHistograms","****   SWAPPING HISTOGRAMS             ****");
1863   Info("SwapGeneratingHistograms","*******************************************");
1864   Info("SwapGeneratingHistograms","*******************************************");
1865   Info("SwapGeneratingHistograms","*******************************************");
1866
1867
1868   TH3D* h = fQSignal;
1869   fQSignal = fQSecondSignal;
1870   fQSecondSignal = h;
1871   fQSecondSignal->Reset();
1872   fQSecondSignal->SetDirectory(0x0);
1873
1874   h = fQBackground;
1875   fQBackground = fQSecondBackground;
1876   fQSecondBackground = h;
1877   fQSecondBackground->Reset();
1878   fQSecondBackground->SetDirectory(0x0);
1879
1880   fSwapped = kTRUE;
1881
1882 }