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