STEER/AliGenInfo.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////
  18 /*
  19
  20 Origin: marian.ivanov@cern.ch
  21
  22 Macro to generate comples MC information - used for Comparison later on
  23 How to use it?
  24
  25 .L $ALICE_ROOT/STEER/AliGenInfo.C+
  26 AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",1,0)
  27 t->Exec();
  28
  29 */
  30
  31 #if !defined(__CINT__) || defined(__MAKECINT__)
  32 #include <stdio.h>
  33 #include <string.h>
  34 //ROOT includes
  35 #include "Rtypes.h"
  36 #include "TFile.h"
  37 #include "TTree.h"
  38 #include "TChain.h"
  39 #include "TCut.h"
  40 #include "TString.h"
  41 #include "TBenchmark.h"
  42 #include "TStopwatch.h"
  43 #include "TParticle.h"
  44 #include "TSystem.h"
  45 #include "TTimer.h"
  46 #include "TVector3.h"
  47 #include "TH1F.h"
  48 #include "TH2F.h"
  49 #include "TCanvas.h"
  50 #include "TPad.h"
  51 #include "TF1.h"
  52 #include "TView.h"
  53 #include "TGeometry.h"
  54 #include "TPolyLine3D.h"
  55
  56 //ALIROOT includes
  57 #include "AliRun.h"
  58 #include "AliStack.h"
  59 #include "AliSimDigits.h"
  60 #include "AliTPCParam.h"
  61 #include "AliTPC.h"
  62 #include "AliTPCLoader.h"
  63 #include "AliDetector.h"
  64 #include "AliTrackReference.h"
  65 #include "AliTPCParamSR.h"
  66 #include "AliTracker.h"
  67 #include "AliMagF.h"
  68 #include "AliHelix.h"
  69 #include "AliPoints.h"
  70
  71 #endif
  72 #include "AliGenInfo.h"
  73 //
  74 //
  75
  76 AliTPCParam * GetTPCParam(){
  77   AliTPCParamSR * par = new AliTPCParamSR;
  78   par->Update();
  79   return par;
  80 }
  81
  82
  83 //_____________________________________________________________________________
  84 Float_t TPCBetheBloch(Float_t bg)
  85 {
  86   //
  87   // Bethe-Bloch energy loss formula
  88   //
  89   const Double_t kp1=0.76176e-1;
  90   const Double_t kp2=10.632;
  91   const Double_t kp3=0.13279e-4;
  92   const Double_t kp4=1.8631;
  93   const Double_t kp5=1.9479;
  94
  95   Double_t dbg = (Double_t) bg;
  96
  97   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
  98
  99   Double_t aa = TMath::Power(beta,kp4);
 100   Double_t bb = TMath::Power(1./dbg,kp5);
 101
 102   bb=TMath::Log(kp3+bb);
 103
 104   return ((Float_t)((kp2-aa-bb)*kp1/aa));
 105 }
 106
 107 AliPointsMI::AliPointsMI(){
 108   fN=0;
 109   fX=0;
 110   fY=0;
 111   fZ=0;
 112   fCapacity = 0;
 113   fLabel0=0;
 114   fLabel1=0;
 115 }
 116
 117
 118 AliPointsMI::AliPointsMI(Int_t n, Float_t *x,Float_t *y, Float_t *z){
 119   //
 120   //
 121   fN=n;
 122   fCapacity = 1000+2*n;
 123   fX= new Float_t[n];
 124   fY= new Float_t[n];
 125   fZ= new Float_t[n];
 126   memcpy(fX,x,n*sizeof(Float_t));
 127   memcpy(fY,y,n*sizeof(Float_t));
 128   memcpy(fZ,z,n*sizeof(Float_t));
 129   fLabel0=0;
 130   fLabel1=0;
 131 }
 132
 133 void AliPointsMI::Reset()
 134 {
 135   fN=0;
 136 }
 137
 138 void AliPointsMI::Reset(AliDetector * det, Int_t particle){
 139   //
 140   // write points from detector points
 141   //
 142   Reset();
 143   if (!det) return;
 144   TObjArray *array = det->Points();
 145   if (!array) return;
 146   for (Int_t i=0;i<array->GetEntriesFast();i++){
 147     AliPoints * points = (AliPoints*) array->At(i);
 148     if (!points) continue;
 149     if (points->GetIndex()!=particle) continue;
 150     Int_t npoints = points->GetN();
 151     if (npoints<2) continue;
 152     Int_t delta = npoints/100;
 153     if (delta<1) delta=1;
 154     if (delta>10) delta=10;
 155     Int_t mypoints = npoints/delta;
 156     //
 157     fN = mypoints;
 158     if (fN>fCapacity){
 159       fCapacity = 1000+2*fN;
 160       delete []fX;
 161       delete []fY;
 162       delete []fZ;
 163       fX = new Float_t[fCapacity];
 164       fY = new Float_t[fCapacity];
 165       fZ = new Float_t[fCapacity];
 166     }
 167     Float_t *xyz = points->GetP();
 168     for (Int_t ipoint=0;ipoint<mypoints;ipoint++){
 169       Int_t index = 3*ipoint*delta;
 170       fX[ipoint]=0;
 171       fY[ipoint]=0;
 172       fZ[ipoint]=0;
 173       if (index+2<npoints*3){
 174         fX[ipoint] = xyz[index];
 175         fY[ipoint] = xyz[index+1];
 176         fZ[ipoint] = xyz[index+2];
 177       }
 178     }
 179   }
 180   fLabel0 = particle;
 181 }
 182
 183
 184 AliPointsMI::~AliPointsMI(){
 185   fN=0;
 186   fCapacity =0;
 187   delete[] fX;
 188   delete[]fY;
 189   delete []fZ;
 190   fX=0;fY=0;fZ=0;
 191 }
 192
 193
 194
 195
 196
 197 ////////////////////////////////////////////////////////////////////////
 198 AliMCInfo::AliMCInfo()
 199 {
 200   fTPCReferences  = new TClonesArray("AliTrackReference",10);
 201   fITSReferences  = new TClonesArray("AliTrackReference",10);
 202   fTRDReferences  = new TClonesArray("AliTrackReference",10);
 203   fTOFReferences  = new TClonesArray("AliTrackReference",10);
 204   fTRdecay.SetTrack(-1);
 205 }
 206
 207 AliMCInfo::~AliMCInfo()
 208 {
 209   if (fTPCReferences) {
 210     delete fTPCReferences;
 211   }
 212   if (fITSReferences){
 213     delete fITSReferences;
 214   }
 215   if (fTRDReferences){
 216     delete fTRDReferences;
 217   }
 218   if (fTOFReferences){
 219     delete fTOFReferences;
 220   }
 221
 222 }
 223
 224
 225
 226 void AliMCInfo::Update()
 227 {
 228   //
 229   //
 230   fMCtracks =1;
 231   if (!fTPCReferences) {
 232     fNTPCRef =0;
 233     return;
 234   }
 235   Float_t direction=1;
 236   //Float_t rlast=0;
 237   fNTPCRef = fTPCReferences->GetEntriesFast();
 238   fNITSRef = fITSReferences->GetEntriesFast();
 239   fNTRDRef = fTRDReferences->GetEntriesFast();
 240   fNTOFRef = fTOFReferences->GetEntriesFast();
 241
 242   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
 243     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
 244     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
 245     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
 246     if (iref==0) direction = newdirection;
 247     if ( newdirection*direction<0){
 248       //changed direction
 249       direction = newdirection;
 250       fMCtracks+=1;
 251     }
 252     //rlast=r;
 253   }
 254   //
 255   // decay info
 256   fTPCdecay=kFALSE;
 257   if (fTRdecay.GetTrack()>0){
 258     fDecayCoord[0] = fTRdecay.X();
 259     fDecayCoord[1] = fTRdecay.Y();
 260     fDecayCoord[2] = fTRdecay.Z();
 261     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
 262       fTPCdecay=kTRUE;
 263     }
 264     else{
 265       fDecayCoord[0] = 0;
 266       fDecayCoord[1] = 0;
 267       fDecayCoord[2] = 0;
 268     }
 269   }
 270   //
 271   //
 272   //digits information update
 273   fRowsWithDigits    = fTPCRow.RowsOn();
 274   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
 275   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
 276   //
 277   //
 278   // calculate primary ionization per cm
 279   if (fParticle.GetPDG()){
 280     fMass = fParticle.GetMass();
 281     fCharge = fParticle.GetPDG()->Charge();
 282     if (fTPCReferences->GetEntriesFast()>0){
 283       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
 284     }
 285     if (fMass>0){
 286       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
 287                               fTrackRef.Py()*fTrackRef.Py()+
 288                               fTrackRef.Pz()*fTrackRef.Pz());
 289       if (p>0.001){
 290         Float_t betagama = p /fMass;
 291         fPrim = TPCBetheBloch(betagama);
 292       }else fPrim=0;
 293     }
 294   }else{
 295     fMass =0;
 296     fPrim =0;
 297   }
 298 }
 299
 300 /////////////////////////////////////////////////////////////////////////////////
 301 /*
 302 void AliGenV0Info::Update()
 303 {
 304   fMCPd[0] = fMCd.fParticle.Px();
 305   fMCPd[1] = fMCd.fParticle.Py();
 306   fMCPd[2] = fMCd.fParticle.Pz();
 307   fMCPd[3] = fMCd.fParticle.P();
 308   //
 309   fMCX[0]  = fMCd.fParticle.Vx();
 310   fMCX[1]  = fMCd.fParticle.Vy();
 311   fMCX[2]  = fMCd.fParticle.Vz();
 312   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 313   //
 314   fPdg[0]    = fMCd.fParticle.GetPdgCode();
 315   fPdg[1]    = fMCm.fParticle.GetPdgCode();
 316   //
 317   fLab[0]    = fMCd.fParticle.GetUniqueID();
 318   fLab[1]    = fMCm.fParticle.GetUniqueID();
 319   //
 320 }
 321 */
 322
 323 void AliGenV0Info::Update(Float_t vertex[3])
 324 {
 325   fMCPd[0] = fMCd.fParticle.Px();
 326   fMCPd[1] = fMCd.fParticle.Py();
 327   fMCPd[2] = fMCd.fParticle.Pz();
 328   fMCPd[3] = fMCd.fParticle.P();
 329   //
 330   fMCX[0]  = fMCd.fParticle.Vx();
 331   fMCX[1]  = fMCd.fParticle.Vy();
 332   fMCX[2]  = fMCd.fParticle.Vz();
 333   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 334   //
 335   fPdg[0]    = fMCd.fParticle.GetPdgCode();
 336   fPdg[1]    = fMCm.fParticle.GetPdgCode();
 337   //
 338   fLab[0]    = fMCd.fParticle.GetUniqueID();
 339   fLab[1]    = fMCm.fParticle.GetUniqueID();
 340   //
 341   //
 342   //
 343   Double_t x1[3],p1[3];
 344   TParticle & pdaughter = fMCd.fParticle;
 345   x1[0] = pdaughter.Vx();
 346   x1[1] = pdaughter.Vy();
 347   x1[2] = pdaughter.Vz();
 348   p1[0] = pdaughter.Px();
 349   p1[1] = pdaughter.Py();
 350   p1[2] = pdaughter.Pz();
 351   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
 352   AliHelix dhelix1(x1,p1,sign);
 353   //
 354   //
 355   Double_t x2[3],p2[3];
 356   //
 357   TParticle & pmother = fMCm.fParticle;
 358   x2[0] = pmother.Vx();
 359   x2[1] = pmother.Vy();
 360   x2[2] = pmother.Vz();
 361   p2[0] = pmother.Px();
 362   p2[1] = pmother.Py();
 363   p2[2] = pmother.Pz();
 364   //
 365   //
 366   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
 367   AliHelix mhelix(x2,p2,sign);
 368
 369   //
 370   //
 371   //
 372   //find intersection linear
 373   //
 374   Double_t distance1, distance2;
 375   Double_t phase[2][2],radius[2];
 376   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 377   Double_t delta1=10000,delta2=10000;
 378   if (points>0){
 379     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 380     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 381     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 382   }
 383   else{
 384     fInvMass=-1;
 385     return;
 386   }
 387   if (points==2){
 388     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 389     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 390     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 391   }
 392   distance1 = TMath::Min(delta1,delta2);
 393   //
 394   //find intersection parabolic
 395   //
 396   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 397   delta1=10000,delta2=10000;
 398
 399   if (points>0){
 400     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 401   }
 402   if (points==2){
 403     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 404   }
 405
 406   distance2 = TMath::Min(delta1,delta2);
 407   //
 408   if (delta1<delta2){
 409     //get V0 info
 410     dhelix1.Evaluate(phase[0][0],fMCXr);
 411     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 412     mhelix.GetMomentum(phase[0][1],fMCPm);
 413     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 414     fMCRr = TMath::Sqrt(radius[0]);
 415   }
 416   else{
 417     dhelix1.Evaluate(phase[1][0],fMCXr);
 418     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 419     mhelix.GetMomentum(phase[1][1],fMCPm);
 420     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 421     fMCRr = TMath::Sqrt(radius[1]);
 422   }
 423   //
 424   //
 425   fMCDist1 = TMath::Sqrt(distance1);
 426   fMCDist2 = TMath::Sqrt(distance2);
 427   //
 428   //
 429   //
 430   Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
 431   //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
 432   Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
 433   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 434   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
 435   vnorm2 = TMath::Sqrt(vnorm2);
 436   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 437   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 438   pnorm2 = TMath::Sqrt(pnorm2);
 439   //
 440   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 441   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 442   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 443   Float_t mass1 = fMCd.fMass;
 444   Float_t mass2 = fMCm.fMass;
 445
 446
 447   Float_t e1    = TMath::Sqrt(mass1*mass1+
 448                               fMCPdr[0]*fMCPdr[0]+
 449                               fMCPdr[1]*fMCPdr[1]+
 450                               fMCPdr[2]*fMCPdr[2]);
 451   Float_t e2    = TMath::Sqrt(mass2*mass2+
 452                               fMCPm[0]*fMCPm[0]+
 453                               fMCPm[1]*fMCPm[1]+
 454                               fMCPm[2]*fMCPm[2]);
 455
 456   fInvMass =
 457     (fMCPm[0]+fMCPdr[0])*(fMCPm[0]+fMCPdr[0])+
 458     (fMCPm[1]+fMCPdr[1])*(fMCPm[1]+fMCPdr[1])+
 459     (fMCPm[2]+fMCPdr[2])*(fMCPm[2]+fMCPdr[2]);
 460
 461   //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
 462   fInvMass = (e1+e2)*(e1+e2)-fInvMass;
 463   if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
 464
 465
 466 }
 467
 468
 469
 470
 471
 472 /////////////////////////////////////////////////////////////////////////////////
 473 void AliGenKinkInfo::Update()
 474 {
 475   fMCPd[0] = fMCd.fParticle.Px();
 476   fMCPd[1] = fMCd.fParticle.Py();
 477   fMCPd[2] = fMCd.fParticle.Pz();
 478   fMCPd[3] = fMCd.fParticle.P();
 479   //
 480   fMCX[0]  = fMCd.fParticle.Vx();
 481   fMCX[1]  = fMCd.fParticle.Vy();
 482   fMCX[2]  = fMCd.fParticle.Vz();
 483   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 484   //
 485   fPdg[0]    = fMCd.fParticle.GetPdgCode();
 486   fPdg[1]    = fMCm.fParticle.GetPdgCode();
 487   //
 488   fLab[0]    = fMCd.fParticle.GetUniqueID();
 489   fLab[1]    = fMCm.fParticle.GetUniqueID();
 490   //
 491   //
 492   //
 493   Double_t x1[3],p1[3];
 494   TParticle & pdaughter = fMCd.fParticle;
 495   x1[0] = pdaughter.Vx();
 496   x1[1] = pdaughter.Vy();
 497   x1[2] = pdaughter.Vz();
 498   p1[0] = pdaughter.Px();
 499   p1[1] = pdaughter.Py();
 500   p1[2] = pdaughter.Pz();
 501   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
 502   AliHelix dhelix1(x1,p1,sign);
 503   //
 504   //
 505   Double_t x2[3],p2[3];
 506   //
 507   TParticle & pmother = fMCm.fParticle;
 508   x2[0] = pmother.Vx();
 509   x2[1] = pmother.Vy();
 510   x2[2] = pmother.Vz();
 511   p2[0] = pmother.Px();
 512   p2[1] = pmother.Py();
 513   p2[2] = pmother.Pz();
 514   //
 515   AliTrackReference & pdecay = fMCm.fTRdecay;
 516   x2[0] = pdecay.X();
 517   x2[1] = pdecay.Y();
 518   x2[2] = pdecay.Z();
 519   p2[0] = pdecay.Px();
 520   p2[1] = pdecay.Py();
 521   p2[2] = pdecay.Pz();
 522   //
 523   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
 524   AliHelix mhelix(x2,p2,sign);
 525
 526   //
 527   //
 528   //
 529   //find intersection linear
 530   //
 531   Double_t distance1, distance2;
 532   Double_t phase[2][2],radius[2];
 533   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 534   Double_t delta1=10000,delta2=10000;
 535   if (points>0){
 536     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 537     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 538     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 539   }
 540   if (points==2){
 541     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 542     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 543     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 544   }
 545   distance1 = TMath::Min(delta1,delta2);
 546   //
 547   //find intersection parabolic
 548   //
 549   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 550   delta1=10000,delta2=10000;
 551
 552   if (points>0){
 553     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 554   }
 555   if (points==2){
 556     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 557   }
 558
 559   distance2 = TMath::Min(delta1,delta2);
 560   //
 561   if (delta1<delta2){
 562     //get V0 info
 563     dhelix1.Evaluate(phase[0][0],fMCXr);
 564     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 565     mhelix.GetMomentum(phase[0][1],fMCPm);
 566     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 567     fMCRr = TMath::Sqrt(radius[0]);
 568   }
 569   else{
 570     dhelix1.Evaluate(phase[1][0],fMCXr);
 571     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 572     mhelix.GetMomentum(phase[1][1],fMCPm);
 573     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 574     fMCRr = TMath::Sqrt(radius[1]);
 575   }
 576   //
 577   //
 578   fMCDist1 = TMath::Sqrt(distance1);
 579   fMCDist2 = TMath::Sqrt(distance2);
 580
 581 }
 582
 583
 584 Float_t AliGenKinkInfo::GetQt(){
 585   //
 586   //
 587   Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
 588   return TMath::Sin(fMCAngle[2])*momentumd;
 589 }
 590
 591
 592
 593
 594 ////////////////////////////////////////////////////////////////////////
 595
 596 ////////////////////////////////////////////////////////////////////////
 597 //
 598 // End of implementation of the class AliMCInfo
 599 //
 600 ////////////////////////////////////////////////////////////////////////
 601
 602
 603
 604 ////////////////////////////////////////////////////////////////////////
 605 digitRow::digitRow()
 606 {
 607   Reset();
 608 }
 609 ////////////////////////////////////////////////////////////////////////
 610 digitRow & digitRow::operator=(const digitRow &digOld)
 611 {
 612   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
 613   return (*this);
 614 }
 615 ////////////////////////////////////////////////////////////////////////
 616 void digitRow::SetRow(Int_t row)
 617 {
 618   if (row >= 8*kgRowBytes) {
 619     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 620     return;
 621   }
 622   Int_t iC = row/8;
 623   Int_t iB = row%8;
 624   SETBIT(fDig[iC],iB);
 625 }
 626
 627 ////////////////////////////////////////////////////////////////////////
 628 Bool_t digitRow::TestRow(Int_t row)
 629 {
 630 //
 631 // return kTRUE if row is on
 632 //
 633   Int_t iC = row/8;
 634   Int_t iB = row%8;
 635   return TESTBIT(fDig[iC],iB);
 636 }
 637 ////////////////////////////////////////////////////////////////////////
 638 Int_t digitRow::RowsOn(Int_t upto)
 639 {
 640 //
 641 // returns number of rows with a digit
 642 // count only rows less equal row number upto
 643 //
 644   Int_t total = 0;
 645   for (Int_t i = 0; i<kgRowBytes; i++) {
 646     for (Int_t j = 0; j < 8; j++) {
 647       if (i*8+j > upto) return total;
 648       if (TESTBIT(fDig[i],j))  total++;
 649     }
 650   }
 651   return total;
 652 }
 653 ////////////////////////////////////////////////////////////////////////
 654 void digitRow::Reset()
 655 {
 656 //
 657 // resets all rows to zero
 658 //
 659   for (Int_t i = 0; i<kgRowBytes; i++) {
 660     fDig[i] <<= 8;
 661   }
 662 }
 663 ////////////////////////////////////////////////////////////////////////
 664 Int_t digitRow::Last()
 665 {
 666 //
 667 // returns the last row number with a digit
 668 // returns -1 if now digits
 669 //
 670   for (Int_t i = kgRowBytes-1; i>=0; i--) {
 671     for (Int_t j = 7; j >= 0; j--) {
 672       if TESTBIT(fDig[i],j) return i*8+j;
 673     }
 674   }
 675   return -1;
 676 }
 677 ////////////////////////////////////////////////////////////////////////
 678 Int_t digitRow::First()
 679 {
 680 //
 681 // returns the first row number with a digit
 682 // returns -1 if now digits
 683 //
 684   for (Int_t i = 0; i<kgRowBytes; i++) {
 685     for (Int_t j = 0; j < 8; j++) {
 686       if (TESTBIT(fDig[i],j)) return i*8+j;
 687     }
 688   }
 689   return -1;
 690 }
 691
 692 ////////////////////////////////////////////////////////////////////////
 693 //
 694 // end of implementation of a class digitRow
 695 //
 696 ////////////////////////////////////////////////////////////////////////
 697
 698 ////////////////////////////////////////////////////////////////////////
 699 AliGenInfoMaker::AliGenInfoMaker()
 700 {
 701   Reset();
 702 }
 703
 704 ////////////////////////////////////////////////////////////////////////
 705 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
 706                                    Int_t nEvents, Int_t firstEvent)
 707 {
 708   Reset();
 709   fFirstEventNr = firstEvent;
 710   fEventNr = firstEvent;
 711   fNEvents = nEvents;
 712   //   fFnRes = fnRes;
 713   sprintf(fFnRes,"%s",fnRes);
 714   //
 715   fLoader = AliRunLoader::Open(fnGalice);
 716   if (gAlice){
 717     delete gAlice->GetRunLoader();
 718     delete gAlice;
 719     gAlice = 0x0;
 720   }
 721   if (fLoader->LoadgAlice()){
 722     cerr<<"Error occured while l"<<endl;
 723   }
 724   Int_t nall = fLoader->GetNumberOfEvents();
 725   if (nEvents==0) {
 726     nEvents =nall;
 727     fNEvents=nall;
 728     fFirstEventNr=0;
 729   }
 730
 731   if (nall<=0){
 732     cerr<<"no events available"<<endl;
 733     fEventNr = 0;
 734     return;
 735   }
 736   if (firstEvent+nEvents>nall) {
 737     fEventNr = nall-firstEvent;
 738     cerr<<"restricted number of events availaible"<<endl;
 739   }
 740   AliMagF * magf = gAlice->Field();
 741   AliTracker::SetFieldMap(magf);
 742 }
 743
 744
 745 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
 746 {
 747   //
 748   if (i<fNParticles) {
 749     if (fGenInfo[i]) return  fGenInfo[i];
 750     fGenInfo[i] = new AliMCInfo;
 751     fNInfos++;
 752     return fGenInfo[i];
 753   }
 754   else
 755     return 0;
 756 }
 757
 758 ////////////////////////////////////////////////////////////////////////
 759 void AliGenInfoMaker::Reset()
 760 {
 761   fEventNr = 0;
 762   fNEvents = 0;
 763   fTreeGenTracks = 0;
 764   fFileGenTracks = 0;
 765   fGenInfo = 0;
 766   fNInfos  = 0;
 767   //
 768   //
 769   fDebug = 0;
 770   fVPrim[0] = -1000.;
 771   fVPrim[1] = -1000.;
 772   fVPrim[2] = -1000.;
 773   fParamTPC = 0;
 774 }
 775 ////////////////////////////////////////////////////////////////////////
 776 AliGenInfoMaker::~AliGenInfoMaker()
 777 {
 778
 779   if (fLoader){
 780     fLoader->UnloadgAlice();
 781     gAlice = 0;
 782     delete fLoader;
 783   }
 784 }
 785
 786 Int_t  AliGenInfoMaker::SetIO()
 787 {
 788   //
 789   //
 790   CreateTreeGenTracks();
 791   if (!fTreeGenTracks) return 1;
 792   //  AliTracker::SetFieldFactor();
 793
 794   fParamTPC = GetTPCParam();
 795   //
 796   return 0;
 797 }
 798
 799 ////////////////////////////////////////////////////////////////////////
 800 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
 801 {
 802   //
 803   //
 804   // SET INPUT
 805   fLoader->SetEventNumber(eventNr);
 806   //
 807   fLoader->LoadHeader();
 808   fLoader->LoadKinematics();
 809   fStack = fLoader->Stack();
 810   //
 811   fLoader->LoadTrackRefs();
 812   fLoader->LoadHits();
 813   fTreeTR = fLoader->TreeTR();
 814   //
 815   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 816   tpcl->LoadDigits();
 817   fTreeD = tpcl->TreeD();
 818   return 0;
 819 }
 820
 821 Int_t AliGenInfoMaker::CloseIOEvent()
 822 {
 823   fLoader->UnloadHeader();
 824   fLoader->UnloadKinematics();
 825   fLoader->UnloadTrackRefs();
 826   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 827   tpcl->UnloadDigits();
 828   return 0;
 829 }
 830
 831 Int_t AliGenInfoMaker::CloseIO()
 832 {
 833   fLoader->UnloadgAlice();
 834   return 0;
 835 }
 836
 837
 838
 839 ////////////////////////////////////////////////////////////////////////
 840 Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
 841 {
 842   fNEvents = nEvents;
 843   fFirstEventNr = firstEventNr;
 844   return Exec();
 845 }
 846
 847 ////////////////////////////////////////////////////////////////////////
 848 Int_t AliGenInfoMaker::Exec()
 849 {
 850   TStopwatch timer;
 851   timer.Start();
 852   Int_t status =SetIO();
 853   if (status>0) return status;
 854   //
 855
 856   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 857        fEventNr++) {
 858     SetIO(fEventNr);
 859     fNParticles = fStack->GetNtrack();
 860     //
 861     fGenInfo = new AliMCInfo*[fNParticles];
 862     for (UInt_t i = 0; i<fNParticles; i++) {
 863       fGenInfo[i]=0;
 864     }
 865     //
 866     cout<<"Start to process event "<<fEventNr<<endl;
 867     cout<<"\tfNParticles = "<<fNParticles<<endl;
 868     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
 869     if (TreeTRLoop()>0) return 1;
 870     //
 871     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
 872     if (TreeDLoop()>0) return 1;
 873     //
 874     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
 875     if (TreeKLoop()>0) return 1;
 876     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 877     //
 878     if (BuildKinkInfo()>0) return 1;
 879     if (BuildV0Info()>0) return 1;
 880     //if (BuildHitLines()>0) return 1;
 881     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 882     //
 883     for (UInt_t i = 0; i<fNParticles; i++) {
 884       if (fGenInfo[i]) delete fGenInfo[i];
 885     }
 886     delete []fGenInfo;
 887     CloseIOEvent();
 888   }
 889   //
 890   CloseIO();
 891   CloseOutputFile();
 892
 893   cerr<<"Exec finished"<<endl;
 894
 895   timer.Stop();
 896   timer.Print();
 897   return 0;
 898 }
 899 ////////////////////////////////////////////////////////////////////////
 900 void AliGenInfoMaker::CreateTreeGenTracks()
 901 {
 902   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 903   if (!fFileGenTracks) {
 904     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 905     return;
 906   }
 907   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 908   AliMCInfo * info = new AliMCInfo;
 909   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
 910   delete info;
 911   //
 912   AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
 913   fTreeKinks = new TTree("genKinksTree","genKinksTree");
 914   fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
 915   delete kinkinfo;
 916   //
 917   AliGenV0Info *v0info = new AliGenV0Info;
 918   fTreeV0 = new TTree("genV0Tree","genV0Tree");
 919   fTreeV0->Branch("MC","AliGenV0Info",&v0info);
 920   delete v0info;
 921   //
 922   //
 923   AliPointsMI * points0 = new AliPointsMI;
 924   AliPointsMI * points1 = new AliPointsMI;
 925   AliPointsMI * points2 = new AliPointsMI;
 926   fTreeHitLines = new TTree("HitLines","HitLines");
 927   fTreeHitLines->Branch("TPC.","AliPointsMI",&points0,32000,10);
 928   fTreeHitLines->Branch("TRD.","AliPointsMI",&points1,32000,10);
 929   fTreeHitLines->Branch("ITS.","AliPointsMI",&points2,32000,10);
 930   Int_t label=0;
 931   fTreeHitLines->Branch("Label",&label,"label/I");
 932   //
 933   fTreeGenTracks->AutoSave();
 934   fTreeKinks->AutoSave();
 935   fTreeV0->AutoSave();
 936   fTreeHitLines->AutoSave();
 937 }
 938 ////////////////////////////////////////////////////////////////////////
 939 void AliGenInfoMaker::CloseOutputFile()
 940 {
 941   if (!fFileGenTracks) {
 942     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 943     return;
 944   }
 945   fFileGenTracks->cd();
 946   fTreeGenTracks->Write();
 947   delete fTreeGenTracks;
 948   fTreeKinks->Write();
 949   delete fTreeKinks;
 950   fTreeV0->Write();
 951   delete fTreeV0;
 952
 953   fFileGenTracks->Close();
 954   delete fFileGenTracks;
 955   return;
 956 }
 957
 958 ////////////////////////////////////////////////////////////////////////
 959 Int_t AliGenInfoMaker::TreeKLoop()
 960 {
 961 //
 962 // open the file with treeK
 963 // loop over all entries there and save information about some tracks
 964 //
 965
 966   AliStack * stack = fStack;
 967   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 968
 969   if (fDebug > 0) {
 970     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 971         <<fEventNr<<endl;
 972   }
 973   Int_t  ipdg = 0;
 974   TParticlePDG *ppdg = 0;
 975   // not all generators give primary vertex position. Take the vertex
 976   // of the particle 0 as primary vertex.
 977   TDatabasePDG  pdg; //get pdg table
 978   //thank you very much root for this
 979   TBranch * br = fTreeGenTracks->GetBranch("MC");
 980   TParticle *particle = stack->ParticleFromTreeK(0);
 981   fVPrim[0] = particle->Vx();
 982   fVPrim[1] = particle->Vy();
 983   fVPrim[2] = particle->Vz();
 984   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 985     // load only particles with TR
 986     AliMCInfo * info = GetInfo(iParticle);
 987     if (!info) continue;
 988     //////////////////////////////////////////////////////////////////////
 989     info->fLabel = iParticle;
 990     //
 991     info->fParticle = *(stack->Particle(iParticle));
 992     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
 993     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
 994     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
 995     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
 996                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
 997     //
 998     //
 999     ipdg = info->fParticle.GetPdgCode();
1000     info->fPdg = ipdg;
1001     ppdg = pdg.GetParticle(ipdg);
1002     info->fEventNr = fEventNr;
1003     info->Update();
1004     //////////////////////////////////////////////////////////////////////
1005     br->SetAddress(&info);
1006     fTreeGenTracks->Fill();
1007   }
1008   fTreeGenTracks->AutoSave();
1009   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1010   return 0;
1011 }
1012
1013
1014 ////////////////////////////////////////////////////////////////////////////////////
1015 Int_t  AliGenInfoMaker::BuildKinkInfo()
1016 {
1017   //
1018   // Build tree with Kink Information
1019   //
1020   AliStack * stack = fStack;
1021   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1022
1023   if (fDebug > 0) {
1024     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1025         <<fEventNr<<endl;
1026   }
1027   //  Int_t  ipdg = 0;
1028   //TParticlePDG *ppdg = 0;
1029   // not all generators give primary vertex position. Take the vertex
1030   // of the particle 0 as primary vertex.
1031   TDatabasePDG  pdg; //get pdg table
1032   //thank you very much root for this
1033   TBranch * br = fTreeKinks->GetBranch("MC");
1034   //
1035   AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
1036   //
1037   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1038     // load only particles with TR
1039     AliMCInfo * info = GetInfo(iParticle);
1040     if (!info) continue;
1041     if (info->fCharge==0) continue;
1042     if (info->fTRdecay.P()<0.13) continue;  //momenta cut
1043     if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
1044     TParticle & particle = info->fParticle;
1045     Int_t first = particle.GetDaughter(0);
1046     if (first ==0) continue;
1047
1048     Int_t last = particle.GetDaughter(1);
1049     if (last ==0) last=first;
1050     AliMCInfo * info2 =  0;
1051     AliMCInfo * dinfo =  0;
1052
1053     for (Int_t id2=first;id2<=last;id2++){
1054       info2 = GetInfo(id2);
1055       if (!info2) continue;
1056       TParticle &p2 = info2->fParticle;
1057       Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
1058       if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
1059       if (!(p2.GetPDG())) continue;
1060       dinfo =info2;
1061
1062       kinkinfo->fMCm = (*info);
1063       kinkinfo->fMCd = (*dinfo);
1064       if (kinkinfo->fMCm.fParticle.GetPDG()==0 ||  kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
1065       kinkinfo->Update();
1066       br->SetAddress(&kinkinfo);
1067       fTreeKinks->Fill();
1068     }
1069     /*
1070     if (dinfo){
1071       kinkinfo->fMCm = (*info);
1072       kinkinfo->fMCd = (*dinfo);
1073       kinkinfo->Update();
1074       br->SetAddress(&kinkinfo);
1075       fTreeKinks->Fill();
1076     }
1077     */
1078   }
1079   fTreeGenTracks->AutoSave();
1080   if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
1081   return 0;
1082 }
1083
1084
1085 ////////////////////////////////////////////////////////////////////////////////////
1086 Int_t  AliGenInfoMaker::BuildV0Info()
1087 {
1088   //
1089   // Build tree with V0 Information
1090   //
1091   AliStack * stack = fStack;
1092   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1093
1094   if (fDebug > 0) {
1095     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1096         <<fEventNr<<endl;
1097   }
1098   //  Int_t  ipdg = 0;
1099   //TParticlePDG *ppdg = 0;
1100   // not all generators give primary vertex position. Take the vertex
1101   // of the particle 0 as primary vertex.
1102   TDatabasePDG  pdg; //get pdg table
1103   //thank you very much root for this
1104   TBranch * br = fTreeV0->GetBranch("MC");
1105   //
1106   AliGenV0Info * v0info = new AliGenV0Info;
1107   //
1108   //
1109   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1110     // load only particles with TR
1111     AliMCInfo * info = GetInfo(iParticle);
1112     if (!info) continue;
1113     if (info->fCharge==0) continue;
1114     //
1115     //
1116     TParticle & particle = info->fParticle;
1117     Int_t mother = particle.GetMother(0);
1118     if (mother <=0) continue;
1119     //
1120     TParticle * motherparticle = stack->Particle(mother);
1121     if (!motherparticle) continue;
1122     //
1123     Int_t last = motherparticle->GetDaughter(1);
1124     if (last==(int)iParticle) continue;
1125     AliMCInfo * info2 =  info;
1126     AliMCInfo * dinfo =  GetInfo(last);
1127     if (!dinfo) continue;
1128     if (!dinfo->fParticle.GetPDG()) continue;
1129     if (!info2->fParticle.GetPDG()) continue;
1130     if (dinfo){
1131       v0info->fMCm = (*info);
1132       v0info->fMCd = (*dinfo);
1133       v0info->Update(fVPrim);
1134       br->SetAddress(&v0info);
1135       fTreeV0->Fill();
1136     }
1137   }
1138   fTreeV0->AutoSave();
1139   if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
1140   return 0;
1141 }
1142
1143
1144
1145
1146 ////////////////////////////////////////////////////////////////////////
1147 Int_t AliGenInfoMaker::BuildHitLines()
1148 {
1149
1150 //
1151 // open the file with treeK
1152 // loop over all entries there and save information about some tracks
1153 //
1154
1155   AliStack * stack = fStack;
1156   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1157
1158   if (fDebug > 0) {
1159     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1160         <<fEventNr<<endl;
1161   }
1162   Int_t  ipdg = 0;
1163   // TParticlePDG *ppdg = 0;
1164   // not all generators give primary vertex position. Take the vertex
1165   // of the particle 0 as primary vertex.
1166   TDatabasePDG  pdg; //get pdg table
1167   //thank you very much root for this
1168   AliPointsMI *tpcp = new AliPointsMI;
1169   AliPointsMI *trdp = new AliPointsMI;
1170   AliPointsMI *itsp = new AliPointsMI;
1171   Int_t label =0;
1172   //
1173   TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
1174   TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");
1175   TBranch * brits = fTreeHitLines->GetBranch("ITS.");
1176   TBranch * brlabel = fTreeHitLines->GetBranch("Label");
1177   brlabel->SetAddress(&label);
1178   brtpc->SetAddress(&tpcp);
1179   brtrd->SetAddress(&trdp);
1180   brits->SetAddress(&itsp);
1181   //
1182   AliDetector *dtpc = gAlice->GetDetector("TPC");
1183   AliDetector *dtrd = gAlice->GetDetector("TRD");
1184   AliDetector *dits = gAlice->GetDetector("ITS");
1185
1186   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1187     // load only particles with TR
1188     AliMCInfo * info = GetInfo(iParticle);
1189     if (!info) continue;
1190     Int_t primpart = info->fPrimPart;
1191     ipdg = info->fParticle.GetPdgCode();
1192     label = iParticle;
1193     //
1194     gAlice->ResetHits();
1195     tpcp->Reset();
1196     itsp->Reset();
1197     trdp->Reset();
1198     tpcp->fLabel1 = ipdg;
1199     trdp->fLabel1 = ipdg;
1200     itsp->fLabel1 = ipdg;
1201     if (dtpc->TreeH()->GetEvent(primpart)){
1202       dtpc->LoadPoints(primpart);
1203       tpcp->Reset(dtpc,iParticle);
1204     }
1205     if (dtrd->TreeH()->GetEvent(primpart)){
1206       dtrd->LoadPoints(primpart);
1207       trdp->Reset(dtrd,iParticle);
1208     }
1209     if (dits->TreeH()->GetEvent(primpart)){
1210       dits->LoadPoints(primpart);
1211       itsp->Reset(dits,iParticle);
1212     }
1213     //
1214     fTreeHitLines->Fill();
1215     dtpc->ResetPoints();
1216     dtrd->ResetPoints();
1217     dits->ResetPoints();
1218   }
1219   delete tpcp;
1220   delete trdp;
1221   delete itsp;
1222   fTreeHitLines->AutoSave();
1223   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1224   return 0;
1225 }
1226
1227
1228 ////////////////////////////////////////////////////////////////////////
1229 Int_t AliGenInfoMaker::TreeDLoop()
1230 {
1231   //
1232   // open the file with treeD
1233   // loop over all entries there and save information about some tracks
1234   //
1235
1236   Int_t nInnerSector = fParamTPC->GetNInnerSector();
1237   Int_t rowShift = 0;
1238   Int_t zero=fParamTPC->GetZeroSup()+6;
1239   //
1240   //
1241   AliSimDigits digitsAddress, *digits=&digitsAddress;
1242   fTreeD->GetBranch("Segment")->SetAddress(&digits);
1243
1244   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
1245   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
1246   for (Int_t i=0; i<sectorsByRows; i++) {
1247     if (!fTreeD->GetEvent(i)) continue;
1248     Int_t sec,row;
1249     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
1250     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
1251     // here I expect that upper sectors follow lower sectors
1252     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
1253     //
1254     digits->ExpandTrackBuffer();
1255     digits->First();
1256     do {
1257       Int_t iRow=digits->CurrentRow();
1258       Int_t iColumn=digits->CurrentColumn();
1259       Short_t digitValue = digits->CurrentDigit();
1260       if (digitValue >= zero) {
1261         Int_t label;
1262         for (Int_t j = 0; j<3; j++) {
1263           //      label = digits->GetTrackID(iRow,iColumn,j);
1264           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
1265           if (label >= (Int_t)fNParticles) { //don't label from bakground event
1266             continue;
1267           }
1268           if (label >= 0 && label <= (Int_t)fNParticles) {
1269             if (fDebug > 6 ) {
1270               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
1271                   <<sec<<" "
1272                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
1273                   <<" "<<row<<endl;
1274             }
1275             AliMCInfo * info = GetInfo(label);
1276             if (info){
1277               info->fTPCRow.SetRow(row+rowShift);
1278             }
1279           }
1280         }
1281       }
1282     } while (digits->Next());
1283   }
1284
1285   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
1286   return 0;
1287 }
1288
1289
1290 ////////////////////////////////////////////////////////////////////////
1291 Int_t AliGenInfoMaker::TreeTRLoop()
1292 {
1293   //
1294   // loop over TrackReferences and store the first one for each track
1295   //
1296   TTree * treeTR = fTreeTR;
1297   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
1298   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
1299   //
1300   //
1301   //track references for TPC
1302   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
1303   TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
1304   TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
1305   TClonesArray* TOFArrayTR = new TClonesArray("AliTrackReference");
1306   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
1307   //
1308   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
1309   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
1310   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
1311   if (treeTR->GetBranch("TOF"))    treeTR->GetBranch("TOF")->SetAddress(&TOFArrayTR);
1312   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
1313   //
1314   //
1315   //
1316   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
1317     treeTR->GetEntry(iPrimPart);
1318     //
1319     // Loop over TPC references
1320     //
1321     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
1322       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
1323       //
1324       if (trackRef->TestBit(BIT(2))){
1325         //if decay
1326         if (trackRef->P()<fgTPCPtCut) continue;
1327         Int_t label = trackRef->GetTrack();
1328         AliMCInfo * info = GetInfo(label);
1329         if (!info) info = MakeInfo(label);
1330         info->fTRdecay = *trackRef;
1331       }
1332       //
1333       if (trackRef->P()<fgTPCPtCut) continue;
1334       Int_t label = trackRef->GetTrack();
1335       AliMCInfo * info = GetInfo(label);
1336       if (!info) info = MakeInfo(label);
1337       if (!info) continue;
1338       info->fPrimPart =  iPrimPart;
1339       TClonesArray & arr = *(info->fTPCReferences);
1340       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1341     }
1342     //
1343     // Loop over ITS references
1344     //
1345     for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
1346       AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);
1347       //
1348       //
1349       if (trackRef->P()<fgTPCPtCut) continue;
1350       Int_t label = trackRef->GetTrack();
1351       AliMCInfo * info = GetInfo(label);
1352       if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
1353       if (!info) info = MakeInfo(label);
1354       if (!info) continue;
1355       info->fPrimPart =  iPrimPart;
1356       TClonesArray & arr = *(info->fITSReferences);
1357       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1358     }
1359     //
1360     // Loop over TRD references
1361     //
1362     for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
1363       AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);
1364       //
1365       if (trackRef->P()<fgTPCPtCut) continue;
1366       Int_t label = trackRef->GetTrack();
1367       AliMCInfo * info = GetInfo(label);
1368       if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
1369       if (!info) info = MakeInfo(label);
1370       if (!info) continue;
1371       info->fPrimPart =  iPrimPart;
1372       TClonesArray & arr = *(info->fTRDReferences);
1373       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1374     }
1375     //
1376     // Loop over TOF references
1377     //
1378     for (Int_t iTrackRef = 0; iTrackRef < TOFArrayTR->GetEntriesFast(); iTrackRef++) {
1379       AliTrackReference *trackRef = (AliTrackReference*)TOFArrayTR->At(iTrackRef);
1380       Int_t label = trackRef->GetTrack();
1381       AliMCInfo * info = GetInfo(label);
1382       if (!info){
1383         if (trackRef->Pt()<fgTPCPtCut) continue;
1384         if ( (!info) && trackRef->Pt()<fgTOFPtCut) continue;
1385       }
1386       if (!info) info = MakeInfo(label);
1387       if (!info) continue;
1388       info->fPrimPart =  iPrimPart;
1389       TClonesArray & arr = *(info->fTOFReferences);
1390       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1391     }
1392     //
1393     // get dacay position
1394     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
1395       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
1396       //
1397       Int_t label = trackRef->GetTrack();
1398       AliMCInfo * info = GetInfo(label);
1399       if (!info) continue;
1400       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
1401       //      if (TMath::Abs(trackRef.X());
1402       info->fTRdecay = *trackRef;
1403     }
1404   }
1405   //
1406   TPCArrayTR->Delete();
1407   delete  TPCArrayTR;
1408   TRDArrayTR->Delete();
1409   delete  TRDArrayTR;
1410   TOFArrayTR->Delete();
1411   delete  TOFArrayTR;
1412
1413   ITSArrayTR->Delete();
1414   delete  ITSArrayTR;
1415   RunArrayTR->Delete();
1416   delete  RunArrayTR;
1417   //
1418   return 0;
1419 }
1420
1421 ////////////////////////////////////////////////////////////////////////
1422 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
1423                                     AliTPCParam *paramTPC) {
1424
1425   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1426   Int_t index[4];
1427   paramTPC->Transform0to1(x,index);
1428   paramTPC->Transform1to2(x,index);
1429   return x[0];
1430 }
1431 ////////////////////////////////////////////////////////////////////////
1432
1433
1434
1435 TH1F * AliComparisonDraw::DrawXY(const char * chx, const char *chy, const char* selection,
1436                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy, Int_t nBinsRes)
1437 {
1438   //
1439   Double_t* bins = CreateLogBins(nbins, minx, maxx);
1440   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
1441   char cut[1000];
1442   sprintf(cut,"%s&&%s",selection,quality);
1443   char expression[1000];
1444   sprintf(expression,"%s:%s>>hRes2",chy,chx);
1445   fTree->Draw(expression, cut, "groff");
1446   TH1F* hMean=0;
1447   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1448   AliLabelAxes(hRes, chx, chy);
1449   //
1450   delete hRes2;
1451   delete[] bins;
1452   fRes  = hRes;
1453   fMean = hMean;
1454   return hRes;
1455 }
1456
1457 TH1F * AliComparisonDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
1458                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
1459 {
1460   //
1461   Double_t* bins = CreateLogBins(nbins, minx, maxx);
1462   Int_t nBinsRes = 30;
1463   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
1464   char cut[1000];
1465   sprintf(cut,"%s&&%s",selection,quality);
1466   char expression[1000];
1467   sprintf(expression,"%s:%s>>hRes2",chy,chx);
1468   fTree->Draw(expression, cut, "groff");
1469   TH1F* hMean=0;
1470   TH1F* hRes = CreateResHisto(hRes2, &hMean);
1471   AliLabelAxes(hRes, chx, chy);
1472   //
1473   delete hRes2;
1474   delete[] bins;
1475   fRes  = hRes;
1476   fMean = hMean;
1477   return hRes;
1478 }
1479
1480 ///////////////////////////////////////////////////////////////////////////////////
1481 ///////////////////////////////////////////////////////////////////////////////////
1482 TH1F * AliComparisonDraw::Eff(const char *variable, const char* selection, const char * quality,
1483                               Int_t nbins, Float_t min, Float_t max)
1484 {
1485   //
1486   //
1487   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
1488   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
1489   char inputGen[1000];
1490   sprintf(inputGen,"%s>>hGen", variable);
1491   fTree->Draw(inputGen, selection, "groff");
1492   char selectionRec[256];
1493   sprintf(selectionRec, "%s && %s", selection, quality);
1494   char inputRec[1000];
1495   sprintf(inputRec,"%s>>hRec", variable);
1496   fTree->Draw(inputRec, selectionRec, "groff");
1497   //
1498   TH1F* hEff = CreateEffHisto(hGen, hRec);
1499   AliLabelAxes(hEff, variable, "#epsilon [%]");
1500   fRes = hEff;
1501   delete hRec;
1502   delete hGen;
1503   return hEff;
1504 }
1505
1506
1507
1508 ///////////////////////////////////////////////////////////////////////////////////
1509 ///////////////////////////////////////////////////////////////////////////////////
1510 TH1F * AliComparisonDraw::EffLog(const char *variable, const char* selection, const char * quality,
1511                               Int_t nbins, Float_t min, Float_t max)
1512 {
1513   //
1514   //
1515   Double_t* bins = CreateLogBins(nbins, min, max);
1516   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
1517   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
1518   char inputGen[1000];
1519   sprintf(inputGen,"%s>>hGen", variable);
1520   fTree->Draw(inputGen, selection, "groff");
1521   char selectionRec[256];
1522   sprintf(selectionRec, "%s && %s", selection, quality);
1523   char inputRec[1000];
1524   sprintf(inputRec,"%s>>hRec", variable);
1525   fTree->Draw(inputRec, selectionRec, "groff");
1526   //
1527   TH1F* hEff = CreateEffHisto(hGen, hRec);
1528   AliLabelAxes(hEff, variable, "#epsilon [%]");
1529   fRes = hEff;
1530   delete hRec;
1531   delete hGen;
1532   delete[] bins;
1533   return hEff;
1534 }
1535
1536
1537 ///////////////////////////////////////////////////////////////////////////////////
1538 ///////////////////////////////////////////////////////////////////////////////////
1539
1540 Double_t* AliComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
1541 {
1542   Double_t* bins = new Double_t[nBins+1];
1543   bins[0] = xMin;
1544   Double_t factor = pow(xMax/xMin, 1./nBins);
1545   for (Int_t i = 1; i <= nBins; i++)
1546     bins[i] = factor * bins[i-1];
1547   return bins;
1548 }
1549
1550
1551
1552
1553 void AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
1554 {
1555   //
1556   histo->GetXaxis()->SetTitle(xAxisTitle);
1557   histo->GetYaxis()->SetTitle(yAxisTitle);
1558 }
1559
1560
1561 TH1F* AliComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
1562 {
1563   //
1564   Int_t nBins = hGen->GetNbinsX();
1565   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
1566   hEff->SetTitle("");
1567   hEff->SetStats(kFALSE);
1568   hEff->SetMinimum(0.);
1569   hEff->SetMaximum(110.);
1570   //
1571   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
1572     Double_t nGen = hGen->GetBinContent(iBin);
1573     Double_t nRec = hRec->GetBinContent(iBin);
1574     if (nGen > 0) {
1575       Double_t eff = nRec/nGen;
1576       hEff->SetBinContent(iBin, 100. * eff);
1577       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
1578       //      if (error == 0) error = sqrt(0.1/nGen);
1579       //
1580       if (error == 0) error = 0.0001;
1581       hEff->SetBinError(iBin, 100. * error);
1582     } else {
1583       hEff->SetBinContent(iBin, 100. * 0.5);
1584       hEff->SetBinError(iBin, 100. * 0.5);
1585     }
1586   }
1587   return hEff;
1588 }
1589
1590
1591
1592 TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits,
1593                      Bool_t overflowBinFits)
1594 {
1595   TVirtualPad* currentPad = gPad;
1596   TAxis* axis = hRes2->GetXaxis();
1597   Int_t nBins = axis->GetNbins();
1598   TH1F* hRes, *hMean;
1599   if (axis->GetXbins()->GetSize()){
1600     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
1601     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
1602   }
1603   else{
1604     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
1605     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
1606
1607   }
1608   hRes->SetStats(false);
1609   hRes->SetOption("E");
1610   hRes->SetMinimum(0.);
1611   //
1612   hMean->SetStats(false);
1613   hMean->SetOption("E");
1614
1615   // create the fit function
1616   //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
1617   //   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1618   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
1619
1620   fitFunc->SetLineWidth(2);
1621   fitFunc->SetFillStyle(0);
1622   // create canvas for fits
1623   TCanvas* canBinFits = NULL;
1624   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
1625   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
1626   Int_t ny = (nPads-1) / nx + 1;
1627   if (drawBinFits) {
1628     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
1629     if (canBinFits) delete canBinFits;
1630     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
1631     canBinFits->Divide(nx, ny);
1632   }
1633
1634   // loop over x bins and fit projection
1635   Int_t dBin = ((overflowBinFits) ? 1 : 0);
1636   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
1637     if (drawBinFits) canBinFits->cd(bin + dBin);
1638     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
1639     //
1640     if (hBin->GetEntries() > 5) {
1641       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
1642       hBin->Fit(fitFunc,"s");
1643       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
1644
1645       if (sigma > 0.){
1646         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
1647         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
1648       }
1649       else{
1650         hRes->SetBinContent(bin, 0.);
1651         hMean->SetBinContent(bin,0);
1652       }
1653       hRes->SetBinError(bin, fitFunc->GetParError(2));
1654       hMean->SetBinError(bin, fitFunc->GetParError(1));
1655
1656       //
1657       //
1658
1659     } else {
1660       hRes->SetBinContent(bin, 0.);
1661       hRes->SetBinError(bin, 0.);
1662       hMean->SetBinContent(bin, 0.);
1663       hMean->SetBinError(bin, 0.);
1664     }
1665
1666
1667     if (drawBinFits) {
1668       char name[256];
1669       if (bin == 0) {
1670         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1671       } else if (bin == nBins+1) {
1672         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1673       } else {
1674         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1675                 axis->GetTitle(), axis->GetBinUpEdge(bin));
1676       }
1677       canBinFits->cd(bin + dBin);
1678       hBin->SetTitle(name);
1679       hBin->SetStats(kTRUE);
1680       hBin->DrawCopy("E");
1681       canBinFits->Update();
1682       canBinFits->Modified();
1683       canBinFits->Update();
1684     }
1685
1686     delete hBin;
1687   }
1688
1689   delete fitFunc;
1690   currentPad->cd();
1691   *phMean = hMean;
1692   return hRes;
1693 }
1694
1695
1696 void   AliComparisonDraw::DrawFriend2D(const char * chx, const char *chy, const char* selection, TTree * tfriend)
1697 {
1698
1699 }
1700
1701
1702 void   AliComparisonDraw::GetPoints3D(const char * label, const char * chpoints, const char* selection, TTree * tpoints, Int_t color,Float_t rmin){
1703   //
1704   //
1705   if (!fPoints) fPoints = new TObjArray;
1706   if (tpoints->GetIndex()==0) tpoints->BuildIndex("Label","Label");
1707   AliPointsMI * points = new AliPointsMI;
1708   TBranch * br = tpoints->GetBranch(chpoints);
1709   br->SetAddress(&points);
1710   Int_t npoints = fTree->Draw(label,selection);
1711   Float_t xyz[30000];
1712   rmin*=rmin;
1713   for (Int_t i=0;i<npoints;i++){
1714     Int_t index = (Int_t)fTree->GetV1()[i];
1715     tpoints->GetEntryWithIndex(index,index);
1716     if (points->fN<2) continue;
1717     Int_t goodpoints=0;
1718     for (Int_t i=0;i<points->fN; i++){
1719       xyz[goodpoints*3]   = points->fX[i];
1720       xyz[goodpoints*3+1] = points->fY[i];
1721       xyz[goodpoints*3+2] = points->fZ[i];
1722       if ( points->fX[i]*points->fX[i]+points->fY[i]*points->fY[i]>rmin) goodpoints++;
1723     }
1724     TPolyMarker3D * marker = new TPolyMarker3D(goodpoints,xyz);
1725     //    marker->SeWidth(1);
1726     marker->SetMarkerColor(color);
1727     marker->SetMarkerStyle(1);
1728     fPoints->AddLast(marker);
1729   }
1730 }
1731
1732 void AliComparisonDraw::Draw3D(Int_t min, Int_t max){
1733   if (!fPoints) return;
1734   Int_t npoints = fPoints->GetEntriesFast();
1735   max = TMath::Min(npoints,max);
1736   min = TMath::Min(npoints,min);
1737
1738   for (Int_t i =min; i<max;i++){
1739     TObject * obj = fPoints->At(i);
1740     if (obj) obj->Draw();
1741   }
1742 }
1743
1744 void AliComparisonDraw:: SavePoints(const char* name)
1745 {
1746   char fname[25];
1747   sprintf(fname,"TH_%s.root",name);
1748   TFile f(fname,"new");
1749   fPoints->Write("Tracks",TObject::kSingleKey);
1750   f.Close();
1751   sprintf(fname,"TH%s.gif",name);
1752   fCanvas->SaveAs(fname);
1753 }
1754
1755 void   AliComparisonDraw::InitView(){
1756   //
1757   //
1758   AliRunLoader* rl = AliRunLoader::Open();
1759   rl->GetEvent(0);
1760   rl->CdGAFile();
1761   //
1762   fCanvas = new TCanvas("cdisplay", "Cluster display",0,0,700,730);
1763   fView   = new TView(1);
1764   fView->SetRange(-330,-360,-330, 360,360, 330);
1765   fCanvas->Clear();
1766   fCanvas->SetFillColor(1);
1767   fCanvas->SetTheta(90.);
1768   fCanvas->SetPhi(0.);
1769   //
1770   TGeometry *geom =(TGeometry*)gDirectory->Get("AliceGeom");
1771   geom->Draw("same");
1772
1773 }