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