PWG1/AliGenInfo.cxx

   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 Generate complex MC information - used for Comparison later on
  23 How to use it?
  24
  25 gSystem->Load("libPWG1.so")
  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 "TROOT.h"
  36 #include "Rtypes.h"
  37 #include "TFile.h"
  38 #include "TTree.h"
  39 #include "TChain.h"
  40 #include "TCut.h"
  41 #include "TString.h"
  42 #include "TStopwatch.h"
  43 #include "TParticle.h"
  44 #include "TSystem.h"
  45 #include "TCanvas.h"
  46 #include "TGeometry.h"
  47 #include "TPolyLine3D.h"
  48
  49 //ALIROOT includes
  50 #include "AliRun.h"
  51 #include "AliStack.h"
  52 #include "AliSimDigits.h"
  53 #include "AliTPCParam.h"
  54 #include "AliTPC.h"
  55 #include "AliTPCLoader.h"
  56 #include "AliDetector.h"
  57 #include "AliTrackReference.h"
  58 #include "AliTPCParamSR.h"
  59 #include "AliTracker.h"
  60 #include "AliMagF.h"
  61 #include "AliHelix.h"
  62 #include "AliTrackPointArray.h"
  63
  64 #endif
  65 #include "AliGenInfo.h"
  66 //
  67 //
  68
  69 ClassImp(AliTPCdigitRow);
  70 ClassImp(AliMCInfo);
  71 ClassImp(AliGenV0Info)
  72 ClassImp(AliGenKinkInfo)
  73 ClassImp(AliGenInfoMaker)
  74
  75
  76
  77 AliTPCParam * GetTPCParam(){
  78   AliTPCParamSR * par = new AliTPCParamSR;
  79   par->Update();
  80   return par;
  81 }
  82
  83
  84 //_____________________________________________________________________________
  85 Float_t TPCBetheBloch(Float_t bg)
  86 {
  87   //
  88   // Bethe-Bloch energy loss formula
  89   //
  90   const Double_t kp1=0.76176e-1;
  91   const Double_t kp2=10.632;
  92   const Double_t kp3=0.13279e-4;
  93   const Double_t kp4=1.8631;
  94   const Double_t kp5=1.9479;
  95
  96   Double_t dbg = (Double_t) bg;
  97
  98   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
  99
 100   Double_t aa = TMath::Power(beta,kp4);
 101   Double_t bb = TMath::Power(1./dbg,kp5);
 102
 103   bb=TMath::Log(kp3+bb);
 104
 105   return ((Float_t)((kp2-aa-bb)*kp1/aa));
 106 }
 107
 108
 109
 110
 111
 112 ////////////////////////////////////////////////////////////////////////
 113 AliMCInfo::AliMCInfo():
 114   fTrackRef(),
 115   fTrackRefOut(),
 116   fTRdecay(),
 117   fPrimPart(0),
 118   fParticle(),
 119   fMass(0),
 120   fCharge(0),
 121   fLabel(0),
 122   fEventNr(),
 123   fMCtracks(),
 124   fPdg(0),
 125   fTPCdecay(0),
 126   fRowsWithDigitsInn(0),
 127   fRowsWithDigits(0),
 128   fRowsTrackLength(0),
 129   fPrim(0),
 130   fTPCRow(),
 131   fNTPCRef(0),                    // tpc references counter
 132   fNITSRef(0),                    // ITS references counter
 133   fNTRDRef(0),                    // TRD references counter
 134   fNTOFRef(0),                    // TOF references counter
 135   fTPCReferences(0),
 136   fITSReferences(0),
 137   fTRDReferences(0),
 138   fTOFReferences(0)
 139 {
 140   fTPCReferences  = new TClonesArray("AliTrackReference",10);
 141   fITSReferences  = new TClonesArray("AliTrackReference",10);
 142   fTRDReferences  = new TClonesArray("AliTrackReference",10);
 143   fTOFReferences  = new TClonesArray("AliTrackReference",10);
 144   fTRdecay.SetTrack(-1);
 145   fCharge = 0;
 146 }
 147
 148 AliMCInfo::AliMCInfo(const AliMCInfo& info):
 149   TObject(),
 150   fTrackRef(info.fTrackRef),
 151   fTrackRefOut(info.fTrackRefOut),
 152   fTRdecay(info.fTRdecay),
 153   fPrimPart(info.fPrimPart),
 154   fParticle(info.fParticle),
 155   fMass(info.fMass),
 156   fCharge(info.fCharge),
 157   fLabel(info.fLabel),
 158   fEventNr(info.fEventNr),
 159   fMCtracks(info.fMCtracks),
 160   fPdg(info.fPdg),
 161   fTPCdecay(info.fTPCdecay),
 162   fRowsWithDigitsInn(info.fRowsWithDigitsInn),
 163   fRowsWithDigits(info.fRowsWithDigits),
 164   fRowsTrackLength(info.fRowsTrackLength),
 165   fPrim(info.fPrim),
 166   fTPCRow(info.fTPCRow),
 167   fNTPCRef(info.fNTPCRef),                    // tpc references counter
 168   fNITSRef(info.fNITSRef),                    // ITS references counter
 169   fNTRDRef(info.fNTRDRef),                    // TRD references counter
 170   fNTOFRef(info.fNTOFRef),                    // TOF references counter
 171   fTPCReferences(0),
 172   fITSReferences(0),
 173   fTRDReferences(0),
 174   fTOFReferences(0)
 175 {
 176   fTPCReferences = (TClonesArray*)info.fTPCReferences->Clone();
 177   fITSReferences = (TClonesArray*)info.fITSReferences->Clone();
 178   fTRDReferences = (TClonesArray*)info.fTRDReferences->Clone();
 179   fTOFReferences = (TClonesArray*)info.fTOFReferences->Clone();
 180 }
 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 = TPCBetheBloch(betagama);
 268       }else fPrim=0;
 269     }
 270   }else{
 271     fMass =0;
 272     fPrim =0;
 273   }
 274 }
 275
 276 /////////////////////////////////////////////////////////////////////////////////
 277 AliGenV0Info::AliGenV0Info():
 278   fMCd(),      //info about daughter particle -
 279   fMCm(),      //info about mother particle   - first particle for V0
 280   fMotherP(),   //particle info about mother particle
 281   fMCDist1(0), //info about closest distance according closest MC - linear DCA
 282   fMCDist2(0),    //info about closest distance parabolic DCA
 283   fMCRr(0),       // rec position of the vertex
 284   fMCR(0),        //exact r position of the vertex
 285   fInvMass(0),  //reconstructed invariant mass -
 286   fPointAngleFi(0), //point angle fi
 287   fPointAngleTh(0), //point angle theta
 288   fPointAngle(0)   //point angle full
 289 {
 290   for (Int_t i=0;i<3; i++){
 291    fMCPdr[i]=0;
 292    fMCX[i]=0;
 293    fMCXr[i]=0;
 294    fMCPm[i]=0;
 295    fMCAngle[i]=0;
 296    fMCPd[i]=0;
 297   }
 298   fMCPd[3]=0;
 299   for (Int_t i=0; i<2;i++){
 300     fPdg[i]=0;
 301     fLab[i]=0;
 302   }
 303 }
 304
 305 void AliGenV0Info::Update(Float_t vertex[3])
 306 {
 307   //
 308   // Update information - calculates derived variables
 309   //
 310
 311   fMCPd[0] = fMCd.fParticle.Px();
 312   fMCPd[1] = fMCd.fParticle.Py();
 313   fMCPd[2] = fMCd.fParticle.Pz();
 314   fMCPd[3] = fMCd.fParticle.P();
 315   //
 316   fMCX[0]  = fMCd.fParticle.Vx();
 317   fMCX[1]  = fMCd.fParticle.Vy();
 318   fMCX[2]  = fMCd.fParticle.Vz();
 319   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 320   //
 321   fPdg[0]    = fMCd.fParticle.GetPdgCode();
 322   fPdg[1]    = fMCm.fParticle.GetPdgCode();
 323   //
 324   fLab[0]    = fMCd.fParticle.GetUniqueID();
 325   fLab[1]    = fMCm.fParticle.GetUniqueID();
 326   //
 327   //
 328   //
 329   Double_t x1[3],p1[3];
 330   TParticle & pdaughter = fMCd.fParticle;
 331   x1[0] = pdaughter.Vx();
 332   x1[1] = pdaughter.Vy();
 333   x1[2] = pdaughter.Vz();
 334   p1[0] = pdaughter.Px();
 335   p1[1] = pdaughter.Py();
 336   p1[2] = pdaughter.Pz();
 337   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
 338   AliHelix dhelix1(x1,p1,sign);
 339   //
 340   //
 341   Double_t x2[3],p2[3];
 342   //
 343   TParticle & pmother = fMCm.fParticle;
 344   x2[0] = pmother.Vx();
 345   x2[1] = pmother.Vy();
 346   x2[2] = pmother.Vz();
 347   p2[0] = pmother.Px();
 348   p2[1] = pmother.Py();
 349   p2[2] = pmother.Pz();
 350   //
 351   //
 352   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
 353   AliHelix mhelix(x2,p2,sign);
 354
 355   //
 356   //
 357   //
 358   //find intersection linear
 359   //
 360   Double_t distance1, distance2;
 361   Double_t phase[2][2],radius[2];
 362   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 363   Double_t delta1=10000,delta2=10000;
 364   if (points>0){
 365     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 366     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 367     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 368   }
 369   else{
 370     fInvMass=-1;
 371     return;
 372   }
 373   if (points==2){
 374     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 375     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 376     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 377   }
 378   distance1 = TMath::Min(delta1,delta2);
 379   //
 380   //find intersection parabolic
 381   //
 382   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 383   delta1=10000,delta2=10000;
 384
 385   if (points>0){
 386     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 387   }
 388   if (points==2){
 389     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 390   }
 391
 392   distance2 = TMath::Min(delta1,delta2);
 393   //
 394   if (delta1<delta2){
 395     //get V0 info
 396     dhelix1.Evaluate(phase[0][0],fMCXr);
 397     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 398     mhelix.GetMomentum(phase[0][1],fMCPm);
 399     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 400     fMCRr = TMath::Sqrt(radius[0]);
 401   }
 402   else{
 403     dhelix1.Evaluate(phase[1][0],fMCXr);
 404     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 405     mhelix.GetMomentum(phase[1][1],fMCPm);
 406     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 407     fMCRr = TMath::Sqrt(radius[1]);
 408   }
 409   //
 410   //
 411   fMCDist1 = TMath::Sqrt(distance1);
 412   fMCDist2 = TMath::Sqrt(distance2);
 413   //
 414   //
 415   //
 416   Float_t v[3] = {fMCXr[0]-vertex[0],fMCXr[1]-vertex[1],fMCXr[2]-vertex[2]};
 417   //Float_t v[3] = {fMCXr[0],fMCXr[1],fMCXr[2]};
 418   Float_t p[3] = {fMCPdr[0]+fMCPm[0], fMCPdr[1]+fMCPm[1],fMCPdr[2]+fMCPm[2]};
 419   Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
 420   Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
 421   vnorm2 = TMath::Sqrt(vnorm2);
 422   Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
 423   Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
 424   pnorm2 = TMath::Sqrt(pnorm2);
 425   //
 426   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
 427   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);
 428   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
 429   Double_t mass1 = fMCd.fMass;
 430   Double_t mass2 = fMCm.fMass;
 431
 432
 433   Double_t e1    = TMath::Sqrt(mass1*mass1+
 434                               fMCPd[0]*fMCPd[0]+
 435                               fMCPd[1]*fMCPd[1]+
 436                               fMCPd[2]*fMCPd[2]);
 437   Double_t e2    = TMath::Sqrt(mass2*mass2+
 438                               fMCPm[0]*fMCPm[0]+
 439                               fMCPm[1]*fMCPm[1]+
 440                               fMCPm[2]*fMCPm[2]);
 441
 442   fInvMass =
 443     (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
 444     (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
 445     (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
 446
 447   //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
 448   fInvMass = (e1+e2)*(e1+e2)-fInvMass;
 449   if (fInvMass>0) fInvMass = TMath::Sqrt(fInvMass);
 450 }
 451
 452
 453
 454 /////////////////////////////////////////////////////////////////////////////////
 455
 456 AliGenKinkInfo::AliGenKinkInfo():
 457   fMCd(),       //info about daughter particle - second particle for V0
 458   fMCm(),       //info about mother particle   - first particle for V0
 459   fMCDist1(0),  //info about closest distance according closest MC - linear DCA
 460   fMCDist2(0),  //info about closest distance parabolic DCA
 461   fMCRr(0),     // rec position of the vertex
 462   fMCR(0)       //exact r position of the vertex
 463 {
 464   //
 465   // default constructor
 466   //
 467   for (Int_t i=0;i<3;i++){
 468     fMCPdr[i]=0;
 469     fMCPd[i]=0;
 470     fMCX[i]=0;
 471     fMCPm[i]=0;
 472     fMCAngle[i]=0;
 473   }
 474   for (Int_t i=0; i<2; i++) {
 475     fPdg[i]= 0;
 476     fLab[i]=0;
 477   }
 478 }
 479
 480 void AliGenKinkInfo::Update()
 481 {
 482   //
 483   // Update information
 484   //  some redundancy - faster acces to the values in analysis code
 485   //
 486   fMCPd[0] = fMCd.fParticle.Px();
 487   fMCPd[1] = fMCd.fParticle.Py();
 488   fMCPd[2] = fMCd.fParticle.Pz();
 489   fMCPd[3] = fMCd.fParticle.P();
 490   //
 491   fMCX[0]  = fMCd.fParticle.Vx();
 492   fMCX[1]  = fMCd.fParticle.Vy();
 493   fMCX[2]  = fMCd.fParticle.Vz();
 494   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 495   //
 496   fPdg[0]    = fMCd.fParticle.GetPdgCode();
 497   fPdg[1]    = fMCm.fParticle.GetPdgCode();
 498   //
 499   fLab[0]    = fMCd.fParticle.GetUniqueID();
 500   fLab[1]    = fMCm.fParticle.GetUniqueID();
 501   //
 502   //
 503   //
 504   Double_t x1[3],p1[3];
 505   TParticle & pdaughter = fMCd.fParticle;
 506   x1[0] = pdaughter.Vx();
 507   x1[1] = pdaughter.Vy();
 508   x1[2] = pdaughter.Vz();
 509   p1[0] = pdaughter.Px();
 510   p1[1] = pdaughter.Py();
 511   p1[2] = pdaughter.Pz();
 512   Double_t sign = (pdaughter.GetPDG()->Charge()>0)? -1:1;
 513   AliHelix dhelix1(x1,p1,sign);
 514   //
 515   //
 516   Double_t x2[3],p2[3];
 517   //
 518   TParticle & pmother = fMCm.fParticle;
 519   x2[0] = pmother.Vx();
 520   x2[1] = pmother.Vy();
 521   x2[2] = pmother.Vz();
 522   p2[0] = pmother.Px();
 523   p2[1] = pmother.Py();
 524   p2[2] = pmother.Pz();
 525   //
 526   AliTrackReference & pdecay = fMCm.fTRdecay;
 527   x2[0] = pdecay.X();
 528   x2[1] = pdecay.Y();
 529   x2[2] = pdecay.Z();
 530   p2[0] = pdecay.Px();
 531   p2[1] = pdecay.Py();
 532   p2[2] = pdecay.Pz();
 533   //
 534   sign = (pmother.GetPDG()->Charge()>0) ? -1:1;
 535   AliHelix mhelix(x2,p2,sign);
 536
 537   //
 538   //
 539   //
 540   //find intersection linear
 541   //
 542   Double_t distance1, distance2;
 543   Double_t phase[2][2],radius[2];
 544   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 545   Double_t delta1=10000,delta2=10000;
 546   if (points>0){
 547     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 548     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 549     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 550   }
 551   if (points==2){
 552     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 553     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 554     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 555   }
 556   distance1 = TMath::Min(delta1,delta2);
 557   //
 558   //find intersection parabolic
 559   //
 560   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
 561   delta1=10000,delta2=10000;
 562
 563   if (points>0){
 564     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
 565   }
 566   if (points==2){
 567     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
 568   }
 569
 570   distance2 = TMath::Min(delta1,delta2);
 571   //
 572   if (delta1<delta2){
 573     //get V0 info
 574     dhelix1.Evaluate(phase[0][0],fMCXr);
 575     dhelix1.GetMomentum(phase[0][0],fMCPdr);
 576     mhelix.GetMomentum(phase[0][1],fMCPm);
 577     dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],fMCAngle);
 578     fMCRr = TMath::Sqrt(radius[0]);
 579   }
 580   else{
 581     dhelix1.Evaluate(phase[1][0],fMCXr);
 582     dhelix1.GetMomentum(phase[1][0],fMCPdr);
 583     mhelix.GetMomentum(phase[1][1],fMCPm);
 584     dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],fMCAngle);
 585     fMCRr = TMath::Sqrt(radius[1]);
 586   }
 587   //
 588   //
 589   fMCDist1 = TMath::Sqrt(distance1);
 590   fMCDist2 = TMath::Sqrt(distance2);
 591
 592 }
 593
 594
 595 Float_t AliGenKinkInfo::GetQt(){
 596   //
 597   //
 598   Float_t momentumd = TMath::Sqrt(fMCPd[0]*fMCPd[0]+fMCPd[1]*fMCPd[1]+fMCPd[2]*fMCPd[2]);
 599   return TMath::Sin(fMCAngle[2])*momentumd;
 600 }
 601
 602
 603
 604
 605 ////////////////////////////////////////////////////////////////////////
 606 AliTPCdigitRow::AliTPCdigitRow()
 607 {
 608   Reset();
 609 }
 610 ////////////////////////////////////////////////////////////////////////
 611 AliTPCdigitRow & AliTPCdigitRow::operator=(const AliTPCdigitRow &digOld)
 612 {
 613   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
 614   return (*this);
 615 }
 616 ////////////////////////////////////////////////////////////////////////
 617 void AliTPCdigitRow::SetRow(Int_t row)
 618 {
 619   if (row >= 8*kgRowBytes) {
 620     cerr<<"AliTPCdigitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 621     return;
 622   }
 623   Int_t iC = row/8;
 624   Int_t iB = row%8;
 625   SETBIT(fDig[iC],iB);
 626 }
 627
 628 ////////////////////////////////////////////////////////////////////////
 629 Bool_t AliTPCdigitRow::TestRow(Int_t row) const
 630 {
 631 //
 632 // return kTRUE if row is on
 633 //
 634   Int_t iC = row/8;
 635   Int_t iB = row%8;
 636   return TESTBIT(fDig[iC],iB);
 637 }
 638 ////////////////////////////////////////////////////////////////////////
 639 Int_t AliTPCdigitRow::RowsOn(Int_t upto) const
 640 {
 641 //
 642 // returns number of rows with a digit
 643 // count only rows less equal row number upto
 644 //
 645   Int_t total = 0;
 646   for (Int_t i = 0; i<kgRowBytes; i++) {
 647     for (Int_t j = 0; j < 8; j++) {
 648       if (i*8+j > upto) return total;
 649       if (TESTBIT(fDig[i],j))  total++;
 650     }
 651   }
 652   return total;
 653 }
 654 ////////////////////////////////////////////////////////////////////////
 655 void AliTPCdigitRow::Reset()
 656 {
 657 //
 658 // resets all rows to zero
 659 //
 660   for (Int_t i = 0; i<kgRowBytes; i++) {
 661     fDig[i] <<= 8;
 662   }
 663 }
 664 ////////////////////////////////////////////////////////////////////////
 665 Int_t AliTPCdigitRow::Last() const
 666 {
 667 //
 668 // returns the last row number with a digit
 669 // returns -1 if now digits
 670 //
 671   for (Int_t i = kgRowBytes-1; i>=0; i--) {
 672     for (Int_t j = 7; j >= 0; j--) {
 673       if TESTBIT(fDig[i],j) return i*8+j;
 674     }
 675   }
 676   return -1;
 677 }
 678 ////////////////////////////////////////////////////////////////////////
 679 Int_t AliTPCdigitRow::First() const
 680 {
 681 //
 682 // returns the first row number with a digit
 683 // returns -1 if now digits
 684 //
 685   for (Int_t i = 0; i<kgRowBytes; i++) {
 686     for (Int_t j = 0; j < 8; j++) {
 687       if (TESTBIT(fDig[i],j)) return i*8+j;
 688     }
 689   }
 690   return -1;
 691 }
 692
 693
 694
 695 ////////////////////////////////////////////////////////////////////////
 696 AliGenInfoMaker::AliGenInfoMaker():
 697   fDebug(0),                   //! debug flag
 698   fEventNr(0),                 //! current event number
 699   fLabel(0),                   //! track label
 700   fNEvents(0),                 //! number of events to process
 701   fFirstEventNr(0),            //! first event to process
 702   fNParticles(0),              //! number of particles in TreeK
 703   fTreeGenTracks(0),           //! output tree with generated tracks
 704   fTreeKinks(0),               //!  output tree with Kinks
 705   fTreeV0(0),                  //!  output tree with V0
 706   fTreeHitLines(0),            //! tree with hit lines
 707   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
 708   fLoader(0),                  //! pointer to the run loader
 709   fTreeD(0),                   //! current tree with digits
 710   fTreeTR(0),                  //! current tree with TR
 711   fStack(0),                   //! current stack
 712   fGenInfo(0),                 //! array with pointers to gen info
 713   fNInfos(0),                  //! number of tracks with infos
 714   fParamTPC(0),                //! AliTPCParam
 715   fTPCPtCut(0.03),
 716   fITSPtCut(0.1),
 717   fTRDPtCut(0.1),
 718   fTOFPtCut(0.1)
 719 {
 720 }
 721
 722 ////////////////////////////////////////////////////////////////////////
 723 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
 724                                  Int_t nEvents, Int_t firstEvent):
 725   fDebug(0),                   //! debug flag
 726   fEventNr(0),                 //! current event number
 727   fLabel(0),                   //! track label
 728   fNEvents(0),                 //! number of events to process
 729   fFirstEventNr(0),            //! first event to process
 730   fNParticles(0),              //! number of particles in TreeK
 731   fTreeGenTracks(0),           //! output tree with generated tracks
 732   fTreeKinks(0),               //!  output tree with Kinks
 733   fTreeV0(0),                  //!  output tree with V0
 734   fTreeHitLines(0),            //! tree with hit lines
 735   fFileGenTracks(0),           //! output file with stored fTreeGenTracks
 736   fLoader(0),                  //! pointer to the run loader
 737   fTreeD(0),                   //! current tree with digits
 738   fTreeTR(0),                  //! current tree with TR
 739   fStack(0),                   //! current stack
 740   fGenInfo(0),                 //! array with pointers to gen info
 741   fNInfos(0),                  //! number of tracks with infos
 742   fParamTPC(0),                //! AliTPCParam
 743   fTPCPtCut(0.03),
 744   fITSPtCut(0.1),
 745   fTRDPtCut(0.1),
 746   fTOFPtCut(0.1)
 747
 748 {
 749   //
 750   //
 751   //
 752   fFirstEventNr = firstEvent;
 753   fEventNr = firstEvent;
 754   fNEvents = nEvents;
 755   sprintf(fFnRes,"%s",fnRes);
 756   //
 757   fLoader = AliRunLoader::Open(fnGalice);
 758   if (gAlice){
 759     delete gAlice->GetRunLoader();
 760     delete gAlice;
 761     gAlice = 0x0;
 762   }
 763   if (fLoader->LoadgAlice()){
 764     cerr<<"Error occured while l"<<endl;
 765   }
 766   Int_t nall = fLoader->GetNumberOfEvents();
 767   if (nEvents==0) {
 768     nEvents =nall;
 769     fNEvents=nall;
 770     fFirstEventNr=0;
 771   }
 772
 773   if (nall<=0){
 774     cerr<<"no events available"<<endl;
 775     fEventNr = 0;
 776     return;
 777   }
 778   if (firstEvent+nEvents>nall) {
 779     fEventNr = nall-firstEvent;
 780     cerr<<"restricted number of events availaible"<<endl;
 781   }
 782   AliMagF * magf = gAlice->Field();
 783   AliTracker::SetFieldMap(magf,0);
 784 }
 785
 786
 787 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
 788 {
 789   //
 790   if (i<fNParticles) {
 791     if (fGenInfo[i]) return  fGenInfo[i];
 792     fGenInfo[i] = new AliMCInfo;
 793     fNInfos++;
 794     return fGenInfo[i];
 795   }
 796   else
 797     return 0;
 798 }
 799
 800 ////////////////////////////////////////////////////////////////////////
 801 AliGenInfoMaker::~AliGenInfoMaker()
 802 {
 803
 804   if (fLoader){
 805     fLoader->UnloadgAlice();
 806     gAlice = 0;
 807     delete fLoader;
 808   }
 809 }
 810
 811 Int_t  AliGenInfoMaker::SetIO()
 812 {
 813   //
 814   //
 815   CreateTreeGenTracks();
 816   if (!fTreeGenTracks) return 1;
 817   //  AliTracker::SetFieldFactor();
 818
 819   fParamTPC = GetTPCParam();
 820   //
 821   return 0;
 822 }
 823
 824 ////////////////////////////////////////////////////////////////////////
 825 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
 826 {
 827   //
 828   //
 829   // SET INPUT
 830   fLoader->SetEventNumber(eventNr);
 831   //
 832   fLoader->LoadHeader();
 833   fLoader->LoadKinematics();
 834   fStack = fLoader->Stack();
 835   //
 836   fLoader->LoadTrackRefs();
 837   fLoader->LoadHits();
 838   fTreeTR = fLoader->TreeTR();
 839   //
 840   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 841   tpcl->LoadDigits();
 842   fTreeD = tpcl->TreeD();
 843   return 0;
 844 }
 845
 846 Int_t AliGenInfoMaker::CloseIOEvent()
 847 {
 848   fLoader->UnloadHeader();
 849   fLoader->UnloadKinematics();
 850   fLoader->UnloadTrackRefs();
 851   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 852   tpcl->UnloadDigits();
 853   return 0;
 854 }
 855
 856 Int_t AliGenInfoMaker::CloseIO()
 857 {
 858   fLoader->UnloadgAlice();
 859   return 0;
 860 }
 861
 862
 863
 864 ////////////////////////////////////////////////////////////////////////
 865 Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
 866 {
 867   fNEvents = nEvents;
 868   fFirstEventNr = firstEventNr;
 869   return Exec();
 870 }
 871
 872 ////////////////////////////////////////////////////////////////////////
 873 Int_t AliGenInfoMaker::Exec()
 874 {
 875   TStopwatch timer;
 876   timer.Start();
 877   Int_t status =SetIO();
 878   if (status>0) return status;
 879   //
 880
 881   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 882        fEventNr++) {
 883     SetIO(fEventNr);
 884     fNParticles = fStack->GetNtrack();
 885     //
 886     fGenInfo = new AliMCInfo*[fNParticles];
 887     for (UInt_t i = 0; i<fNParticles; i++) {
 888       fGenInfo[i]=0;
 889     }
 890     //
 891     cout<<"Start to process event "<<fEventNr<<endl;
 892     cout<<"\tfNParticles = "<<fNParticles<<endl;
 893     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
 894     if (TreeTRLoop()>0) return 1;
 895     //
 896     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
 897     if (TreeDLoop()>0) return 1;
 898     //
 899     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
 900     if (TreeKLoop()>0) return 1;
 901     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 902     //
 903     if (BuildKinkInfo()>0) return 1;
 904     if (BuildV0Info()>0) return 1;
 905     //if (BuildHitLines()>0) return 1;
 906     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 907     //
 908     for (UInt_t i = 0; i<fNParticles; i++) {
 909       if (fGenInfo[i]) delete fGenInfo[i];
 910     }
 911     delete []fGenInfo;
 912     CloseIOEvent();
 913   }
 914   //
 915   CloseIO();
 916   CloseOutputFile();
 917
 918   cerr<<"Exec finished"<<endl;
 919
 920   timer.Stop();
 921   timer.Print();
 922   return 0;
 923 }
 924
 925 ////////////////////////////////////////////////////////////////////////
 926 void AliGenInfoMaker::CreateTreeGenTracks()
 927 {
 928   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 929   if (!fFileGenTracks) {
 930     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 931     return;
 932   }
 933   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 934   AliMCInfo * info = new AliMCInfo;
 935   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
 936   delete info;
 937   //
 938   AliGenKinkInfo *kinkinfo = new AliGenKinkInfo;
 939   fTreeKinks = new TTree("genKinksTree","genKinksTree");
 940   fTreeKinks->Branch("MC","AliGenKinkInfo",&kinkinfo);
 941   delete kinkinfo;
 942   //
 943   AliGenV0Info *v0info = new AliGenV0Info;
 944   fTreeV0 = new TTree("genV0Tree","genV0Tree");
 945   fTreeV0->Branch("MC","AliGenV0Info",&v0info);
 946   delete v0info;
 947   //
 948   //
 949   AliTrackPointArray * points0 = new AliTrackPointArray;
 950   AliTrackPointArray * points1 = new AliTrackPointArray;
 951   AliTrackPointArray * points2 = new AliTrackPointArray;
 952   fTreeHitLines = new TTree("HitLines","HitLines");
 953   fTreeHitLines->Branch("TPC.","AliTrackPointArray",&points0,32000,10);
 954   fTreeHitLines->Branch("TRD.","AliTrackPointArray",&points1,32000,10);
 955   fTreeHitLines->Branch("ITS.","AliTrackPointArray",&points2,32000,10);
 956   Int_t label=0;
 957   fTreeHitLines->Branch("Label",&label,"label/I");
 958   //
 959   fTreeGenTracks->AutoSave();
 960   fTreeKinks->AutoSave();
 961   fTreeV0->AutoSave();
 962   fTreeHitLines->AutoSave();
 963 }
 964 ////////////////////////////////////////////////////////////////////////
 965 void AliGenInfoMaker::CloseOutputFile()
 966 {
 967   if (!fFileGenTracks) {
 968     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 969     return;
 970   }
 971   fFileGenTracks->cd();
 972   fTreeGenTracks->Write();
 973   delete fTreeGenTracks;
 974   fTreeKinks->Write();
 975   delete fTreeKinks;
 976   fTreeV0->Write();
 977   delete fTreeV0;
 978
 979   fFileGenTracks->Close();
 980   delete fFileGenTracks;
 981   return;
 982 }
 983
 984 ////////////////////////////////////////////////////////////////////////
 985 Int_t AliGenInfoMaker::TreeKLoop()
 986 {
 987 //
 988 // open the file with treeK
 989 // loop over all entries there and save information about some tracks
 990 //
 991
 992   AliStack * stack = fStack;
 993   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 994
 995   if (fDebug > 0) {
 996     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 997         <<fEventNr<<endl;
 998   }
 999   Int_t  ipdg = 0;
1000   TParticlePDG *ppdg = 0;
1001   // not all generators give primary vertex position. Take the vertex
1002   // of the particle 0 as primary vertex.
1003   TDatabasePDG  pdg; //get pdg table
1004   //thank you very much root for this
1005   TBranch * br = fTreeGenTracks->GetBranch("MC");
1006   TParticle *particle = stack->ParticleFromTreeK(0);
1007   fVPrim[0] = particle->Vx();
1008   fVPrim[1] = particle->Vy();
1009   fVPrim[2] = particle->Vz();
1010   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1011     // load only particles with TR
1012     AliMCInfo * info = GetInfo(iParticle);
1013     if (!info) continue;
1014     //////////////////////////////////////////////////////////////////////
1015     info->fLabel = iParticle;
1016     //
1017     info->fParticle = *(stack->Particle(iParticle));
1018     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
1019     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
1020     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
1021     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
1022                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
1023     //
1024     //
1025     ipdg = info->fParticle.GetPdgCode();
1026     info->fPdg = ipdg;
1027     ppdg = pdg.GetParticle(ipdg);
1028     info->fEventNr = fEventNr;
1029     info->Update();
1030     //////////////////////////////////////////////////////////////////////
1031     br->SetAddress(&info);
1032     fTreeGenTracks->Fill();
1033   }
1034   fTreeGenTracks->AutoSave();
1035   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1036   return 0;
1037 }
1038
1039
1040 ////////////////////////////////////////////////////////////////////////////////////
1041 Int_t  AliGenInfoMaker::BuildKinkInfo()
1042 {
1043   //
1044   // Build tree with Kink Information
1045   //
1046   AliStack * stack = fStack;
1047   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1048
1049   if (fDebug > 0) {
1050     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1051         <<fEventNr<<endl;
1052   }
1053   //  Int_t  ipdg = 0;
1054   //TParticlePDG *ppdg = 0;
1055   // not all generators give primary vertex position. Take the vertex
1056   // of the particle 0 as primary vertex.
1057   TDatabasePDG  pdg; //get pdg table
1058   //thank you very much root for this
1059   TBranch * br = fTreeKinks->GetBranch("MC");
1060   //
1061   AliGenKinkInfo * kinkinfo = new AliGenKinkInfo;
1062   //
1063   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1064     // load only particles with TR
1065     AliMCInfo * info = GetInfo(iParticle);
1066     if (!info) continue;
1067     if (info->fCharge==0) continue;
1068     if (info->fTRdecay.P()<0.13) continue;  //momenta cut
1069     if (info->fTRdecay.R()>500)  continue;  //R cut - decay outside barrel
1070     TParticle & particle = info->fParticle;
1071     Int_t first = particle.GetDaughter(0);
1072     if (first ==0) continue;
1073
1074     Int_t last = particle.GetDaughter(1);
1075     if (last ==0) last=first;
1076     AliMCInfo * info2 =  0;
1077     AliMCInfo * dinfo =  0;
1078
1079     for (Int_t id2=first;id2<=last;id2++){
1080       info2 = GetInfo(id2);
1081       if (!info2) continue;
1082       TParticle &p2 = info2->fParticle;
1083       Double_t r = TMath::Sqrt(p2.Vx()*p2.Vx()+p2.Vy()*p2.Vy());
1084       if ( TMath::Abs(info->fTRdecay.R()-r)>0.01) continue;
1085       if (!(p2.GetPDG())) continue;
1086       dinfo =info2;
1087
1088       kinkinfo->fMCm = (*info);
1089       kinkinfo->fMCd = (*dinfo);
1090       if (kinkinfo->fMCm.fParticle.GetPDG()==0 ||  kinkinfo->fMCd.fParticle.GetPDG()==0) continue;
1091       kinkinfo->Update();
1092       br->SetAddress(&kinkinfo);
1093       fTreeKinks->Fill();
1094     }
1095     /*
1096     if (dinfo){
1097       kinkinfo->fMCm = (*info);
1098       kinkinfo->fMCd = (*dinfo);
1099       kinkinfo->Update();
1100       br->SetAddress(&kinkinfo);
1101       fTreeKinks->Fill();
1102     }
1103     */
1104   }
1105   fTreeGenTracks->AutoSave();
1106   if (fDebug > 2) cerr<<"end of Kink Loop"<<endl;
1107   return 0;
1108 }
1109
1110
1111 ////////////////////////////////////////////////////////////////////////////////////
1112 Int_t  AliGenInfoMaker::BuildV0Info()
1113 {
1114   //
1115   // Build tree with V0 Information
1116   //
1117   AliStack * stack = fStack;
1118   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1119
1120   if (fDebug > 0) {
1121     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1122         <<fEventNr<<endl;
1123   }
1124   //  Int_t  ipdg = 0;
1125   //TParticlePDG *ppdg = 0;
1126   // not all generators give primary vertex position. Take the vertex
1127   // of the particle 0 as primary vertex.
1128   TDatabasePDG  pdg; //get pdg table
1129   //thank you very much root for this
1130   TBranch * br = fTreeV0->GetBranch("MC");
1131   //
1132   AliGenV0Info * v0info = new AliGenV0Info;
1133   //
1134   //
1135   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1136     // load only particles with TR
1137     AliMCInfo * info = GetInfo(iParticle);
1138     if (!info) continue;
1139     if (info->fCharge==0) continue;
1140     //
1141     //
1142     TParticle & particle = info->fParticle;
1143     Int_t mother = particle.GetMother(0);
1144     if (mother <=0) continue;
1145     //
1146     TParticle * motherparticle = stack->Particle(mother);
1147     if (!motherparticle) continue;
1148     //
1149     Int_t last = motherparticle->GetDaughter(1);
1150     if (last==(int)iParticle) continue;
1151     AliMCInfo * info2 =  info;
1152     AliMCInfo * dinfo =  GetInfo(last);
1153     if (!dinfo) continue;
1154     if (!dinfo->fParticle.GetPDG()) continue;
1155     if (!info2->fParticle.GetPDG()) continue;
1156     if (dinfo){
1157       v0info->fMCm = (*info);
1158       v0info->fMCd = (*dinfo);
1159       v0info->fMotherP = (*motherparticle);
1160       v0info->Update(fVPrim);
1161       br->SetAddress(&v0info);
1162       fTreeV0->Fill();
1163     }
1164   }
1165   fTreeV0->AutoSave();
1166   if (fDebug > 2) cerr<<"end of V0 Loop"<<endl;
1167   return 0;
1168 }
1169
1170
1171
1172
1173 ////////////////////////////////////////////////////////////////////////
1174 Int_t AliGenInfoMaker::BuildHitLines()
1175 {
1176
1177 //
1178 // open the file with treeK
1179 // loop over all entries there and save information about some tracks
1180 //
1181
1182   AliStack * stack = fStack;
1183   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
1184
1185   if (fDebug > 0) {
1186     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
1187         <<fEventNr<<endl;
1188   }
1189 //   Int_t  ipdg = 0;
1190 //   // TParticlePDG *ppdg = 0;
1191 //   // not all generators give primary vertex position. Take the vertex
1192 //   // of the particle 0 as primary vertex.
1193 //   TDatabasePDG  pdg; //get pdg table
1194 //   //thank you very much root for this
1195 //   AliTrackPointArray *tpcp = new AliTrackPointArray;
1196 //   AliTrackPointArray *trdp = new AliTrackPointArray;
1197 //   AliTrackPointArray *itsp = new AliTrackPointArray;
1198 //   Int_t label =0;
1199 //   //
1200 //   TBranch * brtpc = fTreeHitLines->GetBranch("TPC.");
1201 //   TBranch * brtrd = fTreeHitLines->GetBranch("TRD.");
1202 //   TBranch * brits = fTreeHitLines->GetBranch("ITS.");
1203 //   TBranch * brlabel = fTreeHitLines->GetBranch("Label");
1204 //   brlabel->SetAddress(&label);
1205 //   brtpc->SetAddress(&tpcp);
1206 //   brtrd->SetAddress(&trdp);
1207 //   brits->SetAddress(&itsp);
1208 //   //
1209 //   AliDetector *dtpc = gAlice->GetDetector("TPC");
1210 //   AliDetector *dtrd = gAlice->GetDetector("TRD");
1211 //   AliDetector *dits = gAlice->GetDetector("ITS");
1212
1213 //   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
1214 //     // load only particles with TR
1215 //     AliMCInfo * info = GetInfo(iParticle);
1216 //     if (!info) continue;
1217 //     Int_t primpart = info->fPrimPart;
1218 //     ipdg = info->fParticle.GetPdgCode();
1219 //     label = iParticle;
1220 //     //
1221 //     gAlice->ResetHits();
1222 //     tpcp->Reset();
1223 //     itsp->Reset();
1224 //     trdp->Reset();
1225 //     tpcp->fLabel1 = ipdg;
1226 //     trdp->fLabel1 = ipdg;
1227 //     itsp->fLabel1 = ipdg;
1228 //     if (dtpc->TreeH()->GetEvent(primpart)){
1229 //       dtpc->LoadPoints(primpart);
1230 //       tpcp->Reset(dtpc,iParticle);
1231 //     }
1232 //     if (dtrd->TreeH()->GetEvent(primpart)){
1233 //       dtrd->LoadPoints(primpart);
1234 //       trdp->Reset(dtrd,iParticle);
1235 //     }
1236 //     if (dits->TreeH()->GetEvent(primpart)){
1237 //       dits->LoadPoints(primpart);
1238 //       itsp->Reset(dits,iParticle);
1239 //     }
1240 //     //
1241 //     fTreeHitLines->Fill();
1242 //     dtpc->ResetPoints();
1243 //     dtrd->ResetPoints();
1244 //     dits->ResetPoints();
1245 //   }
1246 //   delete tpcp;
1247 //   delete trdp;
1248 //   delete itsp;
1249 //   fTreeHitLines->AutoSave();
1250 //   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
1251   return 0;
1252 }
1253
1254
1255 ////////////////////////////////////////////////////////////////////////
1256 Int_t AliGenInfoMaker::TreeDLoop()
1257 {
1258   //
1259   // open the file with treeD
1260   // loop over all entries there and save information about some tracks
1261   //
1262
1263   Int_t nInnerSector = fParamTPC->GetNInnerSector();
1264   Int_t rowShift = 0;
1265   Int_t zero=fParamTPC->GetZeroSup()+6;
1266   //
1267   //
1268   AliSimDigits digitsAddress, *digits=&digitsAddress;
1269   fTreeD->GetBranch("Segment")->SetAddress(&digits);
1270
1271   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
1272   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
1273   for (Int_t i=0; i<sectorsByRows; i++) {
1274     if (!fTreeD->GetEvent(i)) continue;
1275     Int_t sec,row;
1276     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
1277     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
1278     // here I expect that upper sectors follow lower sectors
1279     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
1280     //
1281     digits->ExpandTrackBuffer();
1282     digits->First();
1283     do {
1284       Int_t iRow=digits->CurrentRow();
1285       Int_t iColumn=digits->CurrentColumn();
1286       Short_t digitValue = digits->CurrentDigit();
1287       if (digitValue >= zero) {
1288         Int_t label;
1289         for (Int_t j = 0; j<3; j++) {
1290           //      label = digits->GetTrackID(iRow,iColumn,j);
1291           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
1292           if (label >= (Int_t)fNParticles) { //don't label from bakground event
1293             continue;
1294           }
1295           if (label >= 0 && label <= (Int_t)fNParticles) {
1296             if (fDebug > 6 ) {
1297               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
1298                   <<sec<<" "
1299                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
1300                   <<" "<<row<<endl;
1301             }
1302             AliMCInfo * info = GetInfo(label);
1303             if (info){
1304               info->fTPCRow.SetRow(row+rowShift);
1305             }
1306           }
1307         }
1308       }
1309     } while (digits->Next());
1310   }
1311
1312   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
1313   return 0;
1314 }
1315
1316
1317 ////////////////////////////////////////////////////////////////////////
1318 Int_t AliGenInfoMaker::TreeTRLoop()
1319 {
1320   //
1321   // loop over TrackReferences and store the first one for each track
1322   //
1323   TTree * treeTR = fTreeTR;
1324   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
1325   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
1326   //
1327   //
1328   //track references for TPC
1329   TClonesArray* tpcArrayTR = new TClonesArray("AliTrackReference");
1330   TClonesArray* itsArrayTR = new TClonesArray("AliTrackReference");
1331   TClonesArray* trdArrayTR = new TClonesArray("AliTrackReference");
1332   TClonesArray* tofArrayTR = new TClonesArray("AliTrackReference");
1333   TClonesArray* runArrayTR = new TClonesArray("AliTrackReference");
1334   //
1335   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&tpcArrayTR);
1336   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&itsArrayTR);
1337   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&trdArrayTR);
1338   if (treeTR->GetBranch("TOF"))    treeTR->GetBranch("TOF")->SetAddress(&tofArrayTR);
1339   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&runArrayTR);
1340   //
1341   //
1342   //
1343   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
1344     treeTR->GetEntry(iPrimPart);
1345     //
1346     // Loop over TPC references
1347     //
1348     for (Int_t iTrackRef = 0; iTrackRef < tpcArrayTR->GetEntriesFast(); iTrackRef++) {
1349       AliTrackReference *trackRef = (AliTrackReference*)tpcArrayTR->At(iTrackRef);
1350       //
1351       if (trackRef->TestBit(BIT(2))){
1352         //if decay
1353         if (trackRef->P()<fTPCPtCut) continue;
1354         Int_t label = trackRef->GetTrack();
1355         AliMCInfo * info = GetInfo(label);
1356         if (!info) info = MakeInfo(label);
1357         info->fTRdecay = *trackRef;
1358       }
1359       //
1360       if (trackRef->P()<fTPCPtCut) continue;
1361       Int_t label = trackRef->GetTrack();
1362       AliMCInfo * info = GetInfo(label);
1363       if (!info) info = MakeInfo(label);
1364       if (!info) continue;
1365       info->fPrimPart =  iPrimPart;
1366       TClonesArray & arr = *(info->fTPCReferences);
1367       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1368     }
1369     //
1370     // Loop over ITS references
1371     //
1372     for (Int_t iTrackRef = 0; iTrackRef < itsArrayTR->GetEntriesFast(); iTrackRef++) {
1373       AliTrackReference *trackRef = (AliTrackReference*)itsArrayTR->At(iTrackRef);
1374       //
1375       //
1376       if (trackRef->P()<fTPCPtCut) continue;
1377       Int_t label = trackRef->GetTrack();
1378       AliMCInfo * info = GetInfo(label);
1379       if ( (!info) && trackRef->Pt()<fITSPtCut) continue;
1380       if (!info) info = MakeInfo(label);
1381       if (!info) continue;
1382       info->fPrimPart =  iPrimPart;
1383       TClonesArray & arr = *(info->fITSReferences);
1384       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1385     }
1386     //
1387     // Loop over TRD references
1388     //
1389     for (Int_t iTrackRef = 0; iTrackRef < trdArrayTR->GetEntriesFast(); iTrackRef++) {
1390       AliTrackReference *trackRef = (AliTrackReference*)trdArrayTR->At(iTrackRef);
1391       //
1392       if (trackRef->P()<fTPCPtCut) continue;
1393       Int_t label = trackRef->GetTrack();
1394       AliMCInfo * info = GetInfo(label);
1395       if ( (!info) && trackRef->Pt()<fTRDPtCut) continue;
1396       if (!info) info = MakeInfo(label);
1397       if (!info) continue;
1398       info->fPrimPart =  iPrimPart;
1399       TClonesArray & arr = *(info->fTRDReferences);
1400       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1401     }
1402     //
1403     // Loop over TOF references
1404     //
1405     for (Int_t iTrackRef = 0; iTrackRef < tofArrayTR->GetEntriesFast(); iTrackRef++) {
1406       AliTrackReference *trackRef = (AliTrackReference*)tofArrayTR->At(iTrackRef);
1407       Int_t label = trackRef->GetTrack();
1408       AliMCInfo * info = GetInfo(label);
1409       if (!info){
1410         if (trackRef->Pt()<fTPCPtCut) continue;
1411         if ( (!info) && trackRef->Pt()<fTOFPtCut) continue;
1412       }
1413       if (!info) info = MakeInfo(label);
1414       if (!info) continue;
1415       info->fPrimPart =  iPrimPart;
1416       TClonesArray & arr = *(info->fTOFReferences);
1417       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
1418     }
1419     //
1420     // get dacay position
1421     for (Int_t iTrackRef = 0; iTrackRef < runArrayTR->GetEntriesFast(); iTrackRef++) {
1422       AliTrackReference *trackRef = (AliTrackReference*)runArrayTR->At(iTrackRef);
1423       //
1424       Int_t label = trackRef->GetTrack();
1425       AliMCInfo * info = GetInfo(label);
1426       if (!info) continue;
1427       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
1428       //      if (TMath::Abs(trackRef.X());
1429       info->fTRdecay = *trackRef;
1430     }
1431   }
1432   //
1433   tpcArrayTR->Delete();
1434   delete  tpcArrayTR;
1435   trdArrayTR->Delete();
1436   delete  trdArrayTR;
1437   tofArrayTR->Delete();
1438   delete  tofArrayTR;
1439   itsArrayTR->Delete();
1440   delete  itsArrayTR;
1441   runArrayTR->Delete();
1442   delete  runArrayTR;
1443   //
1444   return 0;
1445 }
1446
1447 ////////////////////////////////////////////////////////////////////////
1448 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
1449                                     AliTPCParam *paramTPC) const {
1450
1451   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
1452   Int_t index[4];
1453   paramTPC->Transform0to1(x,index);
1454   paramTPC->Transform1to2(x,index);
1455   return x[0];
1456 }
1457 ////////////////////////////////////////////////////////////////////////
1458
1459
1460