STEER/AliGenInfo.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 ///////////////////////////////////////////////////////////////////////////
  18 /*
  19
  20 Origin: marian.ivanov@cern.ch
  21
  22 Macro to generate comples MC information - used for Comparison later on
  23 How to use it?
  24
  25 .L $ALICE_ROOT/STEER/AliGenInfo.C+
  26 AliGenInfoMaker *t = new AliGenInfoMaker("galice.root","genTracks.root",1,0)
  27 t->Exec();
  28
  29 */
  30
  31 #if !defined(__CINT__) || defined(__MAKECINT__)
  32 #include <stdio.h>
  33 #include <string.h>
  34 //ROOT includes
  35 #include "Rtypes.h"
  36 #include "TFile.h"
  37 #include "TTree.h"
  38 #include "TChain.h"
  39 #include "TCut.h"
  40 #include "TString.h"
  41 #include "TBenchmark.h"
  42 #include "TStopwatch.h"
  43 #include "TParticle.h"
  44 #include "TSystem.h"
  45 #include "TTimer.h"
  46 #include "TVector3.h"
  47 #include "TH1F.h"
  48 #include "TH2F.h"
  49 #include "TCanvas.h"
  50 #include "TPad.h"
  51 #include "TF1.h"
  52
  53 //ALIROOT includes
  54 #include "AliRun.h"
  55 #include "AliStack.h"
  56 #include "AliSimDigits.h"
  57 #include "AliTPCParam.h"
  58 #include "AliTPC.h"
  59 #include "AliTPCLoader.h"
  60 #include "AliDetector.h"
  61 #include "AliTrackReference.h"
  62 #include "AliTPCParamSR.h"
  63 #include "AliTracker.h"
  64 #include "AliMagF.h"
  65 #endif
  66 #include "AliGenInfo.h"
  67 //
  68 //
  69
  70 AliTPCParam * GetTPCParam(){
  71   AliTPCParamSR * par = new AliTPCParamSR;
  72   par->Update();
  73   return par;
  74 }
  75
  76
  77 //_____________________________________________________________________________
  78 Float_t TPCBetheBloch(Float_t bg)
  79 {
  80   //
  81   // Bethe-Bloch energy loss formula
  82   //
  83   const Double_t kp1=0.76176e-1;
  84   const Double_t kp2=10.632;
  85   const Double_t kp3=0.13279e-4;
  86   const Double_t kp4=1.8631;
  87   const Double_t kp5=1.9479;
  88
  89   Double_t dbg = (Double_t) bg;
  90
  91   Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
  92
  93   Double_t aa = TMath::Power(beta,kp4);
  94   Double_t bb = TMath::Power(1./dbg,kp5);
  95
  96   bb=TMath::Log(kp3+bb);
  97
  98   return ((Float_t)((kp2-aa-bb)*kp1/aa));
  99 }
 100
 101
 102
 103
 104
 105 ////////////////////////////////////////////////////////////////////////
 106 AliMCInfo::AliMCInfo()
 107 {
 108   fTPCReferences  = new TClonesArray("AliTrackReference",10);
 109   fITSReferences  = new TClonesArray("AliTrackReference",10);
 110   fTRDReferences  = new TClonesArray("AliTrackReference",10);
 111   fTRdecay.SetTrack(-1);
 112 }
 113
 114 AliMCInfo::~AliMCInfo()
 115 {
 116   if (fTPCReferences) {
 117     delete fTPCReferences;
 118   }
 119   if (fITSReferences){
 120     delete fITSReferences;
 121   }
 122   if (fTRDReferences){
 123     delete fTRDReferences;
 124   }
 125 }
 126
 127
 128
 129 void AliMCInfo::Update()
 130 {
 131   //
 132   //
 133   fMCtracks =1;
 134   if (!fTPCReferences) {
 135     fNTPCRef =0;
 136     return;
 137   }
 138   Float_t direction=1;
 139   //Float_t rlast=0;
 140   fNTPCRef = fTPCReferences->GetEntriesFast();
 141   fNITSRef = fITSReferences->GetEntriesFast();
 142
 143   for (Int_t iref =0;iref<fTPCReferences->GetEntriesFast();iref++){
 144     AliTrackReference * ref = (AliTrackReference *) fTPCReferences->At(iref);
 145     //Float_t r = (ref->X()*ref->X()+ref->Y()*ref->Y());
 146     Float_t newdirection = ref->X()*ref->Px()+ref->Y()*ref->Py(); //inside or outside
 147     if (iref==0) direction = newdirection;
 148     if ( newdirection*direction<0){
 149       //changed direction
 150       direction = newdirection;
 151       fMCtracks+=1;
 152     }
 153     //rlast=r;
 154   }
 155   //
 156   // decay info
 157   fTPCdecay=kFALSE;
 158   if (fTRdecay.GetTrack()>0){
 159     fDecayCoord[0] = fTRdecay.X();
 160     fDecayCoord[1] = fTRdecay.Y();
 161     fDecayCoord[2] = fTRdecay.Z();
 162     if ( (fTRdecay.R()<250)&&(fTRdecay.R()>85) && (TMath::Abs(fTRdecay.Z())<250) ){
 163       fTPCdecay=kTRUE;
 164     }
 165     else{
 166       fDecayCoord[0] = 0;
 167       fDecayCoord[1] = 0;
 168       fDecayCoord[2] = 0;
 169     }
 170   }
 171   //
 172   //
 173   //digits information update
 174   fRowsWithDigits    = fTPCRow.RowsOn();
 175   fRowsWithDigitsInn = fTPCRow.RowsOn(63); // 63 = number of inner rows
 176   fRowsTrackLength   = fTPCRow.Last() - fTPCRow.First();
 177   //
 178   //
 179   // calculate primary ionization per cm
 180   if (fParticle.GetPDG()){
 181     fMass = fParticle.GetMass();
 182     fCharge = fParticle.GetPDG()->Charge();
 183     if (fTPCReferences->GetEntriesFast()>0){
 184       fTrackRef = *((AliTrackReference*)fTPCReferences->At(0));
 185     }
 186     if (fMass>0){
 187       Float_t p = TMath::Sqrt(fTrackRef.Px()*fTrackRef.Px()+
 188                               fTrackRef.Py()*fTrackRef.Py()+
 189                               fTrackRef.Pz()*fTrackRef.Pz());
 190       if (p>0.001){
 191         Float_t betagama = p /fMass;
 192         fPrim = TPCBetheBloch(betagama);
 193       }else fPrim=0;
 194     }
 195   }else{
 196     fMass =0;
 197     fPrim =0;
 198   }
 199 }
 200
 201 /////////////////////////////////////////////////////////////////////////////////
 202 void AliGenV0Info::Update()
 203 {
 204   fMCPd[0] = fMCd.fParticle.Px();
 205   fMCPd[1] = fMCd.fParticle.Py();
 206   fMCPd[2] = fMCd.fParticle.Pz();
 207   fMCPd[3] = fMCd.fParticle.P();
 208   fMCX[0]  = fMCd.fParticle.Vx();
 209   fMCX[1]  = fMCd.fParticle.Vy();
 210   fMCX[2]  = fMCd.fParticle.Vz();
 211   fMCR       = TMath::Sqrt( fMCX[0]*fMCX[0]+fMCX[1]*fMCX[1]);
 212   fPdg[0]    = fMCd.fParticle.GetPdgCode();
 213   fPdg[1]    = fMCm.fParticle.GetPdgCode();
 214   //
 215   fLab[0]    = fMCd.fParticle.GetUniqueID();
 216   fLab[1]    = fMCm.fParticle.GetUniqueID();
 217
 218 }
 219
 220
 221
 222
 223
 224 ////////////////////////////////////////////////////////////////////////
 225
 226 ////////////////////////////////////////////////////////////////////////
 227 //
 228 // End of implementation of the class AliMCInfo
 229 //
 230 ////////////////////////////////////////////////////////////////////////
 231
 232
 233
 234 ////////////////////////////////////////////////////////////////////////
 235 digitRow::digitRow()
 236 {
 237   Reset();
 238 }
 239 ////////////////////////////////////////////////////////////////////////
 240 digitRow & digitRow::operator=(const digitRow &digOld)
 241 {
 242   for (Int_t i = 0; i<kgRowBytes; i++) fDig[i] = digOld.fDig[i];
 243   return (*this);
 244 }
 245 ////////////////////////////////////////////////////////////////////////
 246 void digitRow::SetRow(Int_t row)
 247 {
 248   if (row >= 8*kgRowBytes) {
 249     cerr<<"digitRow::SetRow: index "<<row<<" out of bounds."<<endl;
 250     return;
 251   }
 252   Int_t iC = row/8;
 253   Int_t iB = row%8;
 254   SETBIT(fDig[iC],iB);
 255 }
 256
 257 ////////////////////////////////////////////////////////////////////////
 258 Bool_t digitRow::TestRow(Int_t row)
 259 {
 260 //
 261 // return kTRUE if row is on
 262 //
 263   Int_t iC = row/8;
 264   Int_t iB = row%8;
 265   return TESTBIT(fDig[iC],iB);
 266 }
 267 ////////////////////////////////////////////////////////////////////////
 268 Int_t digitRow::RowsOn(Int_t upto)
 269 {
 270 //
 271 // returns number of rows with a digit
 272 // count only rows less equal row number upto
 273 //
 274   Int_t total = 0;
 275   for (Int_t i = 0; i<kgRowBytes; i++) {
 276     for (Int_t j = 0; j < 8; j++) {
 277       if (i*8+j > upto) return total;
 278       if (TESTBIT(fDig[i],j))  total++;
 279     }
 280   }
 281   return total;
 282 }
 283 ////////////////////////////////////////////////////////////////////////
 284 void digitRow::Reset()
 285 {
 286 //
 287 // resets all rows to zero
 288 //
 289   for (Int_t i = 0; i<kgRowBytes; i++) {
 290     fDig[i] <<= 8;
 291   }
 292 }
 293 ////////////////////////////////////////////////////////////////////////
 294 Int_t digitRow::Last()
 295 {
 296 //
 297 // returns the last row number with a digit
 298 // returns -1 if now digits
 299 //
 300   for (Int_t i = kgRowBytes-1; i>=0; i--) {
 301     for (Int_t j = 7; j >= 0; j--) {
 302       if TESTBIT(fDig[i],j) return i*8+j;
 303     }
 304   }
 305   return -1;
 306 }
 307 ////////////////////////////////////////////////////////////////////////
 308 Int_t digitRow::First()
 309 {
 310 //
 311 // returns the first row number with a digit
 312 // returns -1 if now digits
 313 //
 314   for (Int_t i = 0; i<kgRowBytes; i++) {
 315     for (Int_t j = 0; j < 8; j++) {
 316       if (TESTBIT(fDig[i],j)) return i*8+j;
 317     }
 318   }
 319   return -1;
 320 }
 321
 322 ////////////////////////////////////////////////////////////////////////
 323 //
 324 // end of implementation of a class digitRow
 325 //
 326 ////////////////////////////////////////////////////////////////////////
 327
 328 ////////////////////////////////////////////////////////////////////////
 329 AliGenInfoMaker::AliGenInfoMaker()
 330 {
 331   Reset();
 332 }
 333
 334 ////////////////////////////////////////////////////////////////////////
 335 AliGenInfoMaker::AliGenInfoMaker(const char * fnGalice, const char* fnRes,
 336                                    Int_t nEvents, Int_t firstEvent)
 337 {
 338   Reset();
 339   fFirstEventNr = firstEvent;
 340   fEventNr = firstEvent;
 341   fNEvents = nEvents;
 342   //   fFnRes = fnRes;
 343   sprintf(fFnRes,"%s",fnRes);
 344   //
 345   fLoader = AliRunLoader::Open(fnGalice);
 346   if (gAlice){
 347     delete gAlice->GetRunLoader();
 348     delete gAlice;
 349     gAlice = 0x0;
 350   }
 351   if (fLoader->LoadgAlice()){
 352     cerr<<"Error occured while l"<<endl;
 353   }
 354   Int_t nall = fLoader->GetNumberOfEvents();
 355   if (nEvents==0) {
 356     nEvents =nall;
 357     fNEvents=nall;
 358     fFirstEventNr=0;
 359   }
 360
 361   if (nall<=0){
 362     cerr<<"no events available"<<endl;
 363     fEventNr = 0;
 364     return;
 365   }
 366   if (firstEvent+nEvents>nall) {
 367     fEventNr = nall-firstEvent;
 368     cerr<<"restricted number of events availaible"<<endl;
 369   }
 370   AliMagF * magf = gAlice->Field();
 371   AliTracker::SetFieldMap(magf);
 372 }
 373
 374
 375 AliMCInfo * AliGenInfoMaker::MakeInfo(UInt_t i)
 376 {
 377   //
 378   if (i<fNParticles) {
 379     if (fGenInfo[i]) return  fGenInfo[i];
 380     fGenInfo[i] = new AliMCInfo;
 381     fNInfos++;
 382     return fGenInfo[i];
 383   }
 384   else
 385     return 0;
 386 }
 387
 388 ////////////////////////////////////////////////////////////////////////
 389 void AliGenInfoMaker::Reset()
 390 {
 391   fEventNr = 0;
 392   fNEvents = 0;
 393   fTreeGenTracks = 0;
 394   fFileGenTracks = 0;
 395   fGenInfo = 0;
 396   fNInfos  = 0;
 397   //
 398   //
 399   fDebug = 0;
 400   fVPrim[0] = -1000.;
 401   fVPrim[1] = -1000.;
 402   fVPrim[2] = -1000.;
 403   fParamTPC = 0;
 404 }
 405 ////////////////////////////////////////////////////////////////////////
 406 AliGenInfoMaker::~AliGenInfoMaker()
 407 {
 408
 409   if (fLoader){
 410     fLoader->UnloadgAlice();
 411     gAlice = 0;
 412     delete fLoader;
 413   }
 414 }
 415
 416 Int_t  AliGenInfoMaker::SetIO()
 417 {
 418   //
 419   //
 420   CreateTreeGenTracks();
 421   if (!fTreeGenTracks) return 1;
 422   //  AliTracker::SetFieldFactor();
 423
 424   fParamTPC = GetTPCParam();
 425   //
 426   return 0;
 427 }
 428
 429 ////////////////////////////////////////////////////////////////////////
 430 Int_t AliGenInfoMaker::SetIO(Int_t eventNr)
 431 {
 432   //
 433   //
 434   // SET INPUT
 435   fLoader->SetEventNumber(eventNr);
 436   //
 437   fLoader->LoadHeader();
 438   fLoader->LoadKinematics();
 439   fStack = fLoader->Stack();
 440   //
 441   fLoader->LoadTrackRefs();
 442   fTreeTR = fLoader->TreeTR();
 443   //
 444   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 445   tpcl->LoadDigits();
 446   fTreeD = tpcl->TreeD();
 447   return 0;
 448 }
 449
 450 Int_t AliGenInfoMaker::CloseIOEvent()
 451 {
 452   fLoader->UnloadHeader();
 453   fLoader->UnloadKinematics();
 454   fLoader->UnloadTrackRefs();
 455   AliTPCLoader * tpcl = (AliTPCLoader*)fLoader->GetLoader("TPCLoader");
 456   tpcl->UnloadDigits();
 457   return 0;
 458 }
 459
 460 Int_t AliGenInfoMaker::CloseIO()
 461 {
 462   fLoader->UnloadgAlice();
 463   return 0;
 464 }
 465
 466
 467
 468 ////////////////////////////////////////////////////////////////////////
 469 Int_t AliGenInfoMaker::Exec(Int_t nEvents, Int_t firstEventNr)
 470 {
 471   fNEvents = nEvents;
 472   fFirstEventNr = firstEventNr;
 473   return Exec();
 474 }
 475
 476 ////////////////////////////////////////////////////////////////////////
 477 Int_t AliGenInfoMaker::Exec()
 478 {
 479   TStopwatch timer;
 480   timer.Start();
 481   Int_t status =SetIO();
 482   if (status>0) return status;
 483   //
 484
 485   for (fEventNr = fFirstEventNr; fEventNr < fFirstEventNr+fNEvents;
 486        fEventNr++) {
 487     SetIO(fEventNr);
 488     fNParticles = fStack->GetNtrack();
 489     //
 490     fGenInfo = new AliMCInfo*[fNParticles];
 491     for (UInt_t i = 0; i<fNParticles; i++) {
 492       fGenInfo[i]=0;
 493     }
 494     //
 495     cout<<"Start to process event "<<fEventNr<<endl;
 496     cout<<"\tfNParticles = "<<fNParticles<<endl;
 497     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeTR"<<endl;
 498     if (TreeTRLoop()>0) return 1;
 499     //
 500     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeD"<<endl;
 501     if (TreeDLoop()>0) return 1;
 502     //
 503     if (fDebug>2) cout<<"\n\n\n\tStart loop over TreeK"<<endl;
 504     if (TreeKLoop()>0) return 1;
 505     if (fDebug>2) cout<<"\tEnd loop over TreeK"<<endl;
 506     for (UInt_t i = 0; i<fNParticles; i++) {
 507       if (fGenInfo[i]) delete fGenInfo[i];
 508     }
 509     delete []fGenInfo;
 510     CloseIOEvent();
 511   }
 512   //
 513   CloseIO();
 514   CloseOutputFile();
 515
 516   cerr<<"Exec finished"<<endl;
 517
 518   timer.Stop();
 519   timer.Print();
 520   return 0;
 521 }
 522 ////////////////////////////////////////////////////////////////////////
 523 void AliGenInfoMaker::CreateTreeGenTracks()
 524 {
 525   fFileGenTracks = TFile::Open(fFnRes,"RECREATE");
 526   if (!fFileGenTracks) {
 527     cerr<<"Error in CreateTreeGenTracks: cannot open file "<<fFnRes<<endl;
 528     return;
 529   }
 530   fTreeGenTracks = new TTree("genTracksTree","genTracksTree");
 531   AliMCInfo * info = new AliMCInfo;
 532   //
 533   fTreeGenTracks->Branch("MC","AliMCInfo",&info);
 534   delete info;
 535   fTreeGenTracks->AutoSave();
 536 }
 537 ////////////////////////////////////////////////////////////////////////
 538 void AliGenInfoMaker::CloseOutputFile()
 539 {
 540   if (!fFileGenTracks) {
 541     cerr<<"File "<<fFnRes<<" not found as an open file."<<endl;
 542     return;
 543   }
 544   fFileGenTracks->cd();
 545   fTreeGenTracks->Write();
 546   delete fTreeGenTracks;
 547   fFileGenTracks->Close();
 548   delete fFileGenTracks;
 549   return;
 550 }
 551
 552 ////////////////////////////////////////////////////////////////////////
 553 Int_t AliGenInfoMaker::TreeKLoop()
 554 {
 555 //
 556 // open the file with treeK
 557 // loop over all entries there and save information about some tracks
 558 //
 559
 560   AliStack * stack = fStack;
 561   if (!stack) {cerr<<"Stack was not found!\n"; return 1;}
 562
 563   if (fDebug > 0) {
 564     cout<<"There are "<<fNParticles<<" primary and secondary particles in event "
 565         <<fEventNr<<endl;
 566   }
 567   Int_t  ipdg = 0;
 568   TParticlePDG *ppdg = 0;
 569   // not all generators give primary vertex position. Take the vertex
 570   // of the particle 0 as primary vertex.
 571   TDatabasePDG  pdg; //get pdg table
 572   //thank you very much root for this
 573   TBranch * br = fTreeGenTracks->GetBranch("MC");
 574   TParticle *particle = stack->ParticleFromTreeK(0);
 575   fVPrim[0] = particle->Vx();
 576   fVPrim[1] = particle->Vy();
 577   fVPrim[2] = particle->Vz();
 578   for (UInt_t iParticle = 0; iParticle < fNParticles; iParticle++) {
 579     // load only particles with TR
 580     AliMCInfo * info = GetInfo(iParticle);
 581     if (!info) continue;
 582     //////////////////////////////////////////////////////////////////////
 583     info->fLabel = iParticle;
 584     //
 585     info->fParticle = *(stack->Particle(iParticle));
 586     info->fVDist[0] = info->fParticle.Vx()-fVPrim[0];
 587     info->fVDist[1] = info->fParticle.Vy()-fVPrim[1];
 588     info->fVDist[2] = info->fParticle.Vz()-fVPrim[2];
 589     info->fVDist[3] = TMath::Sqrt(info->fVDist[0]*info->fVDist[0]+
 590                                   info->fVDist[1]*info->fVDist[1]+info->fVDist[2]*info->fVDist[2]);
 591     //
 592     //
 593     ipdg = info->fParticle.GetPdgCode();
 594     info->fPdg = ipdg;
 595     ppdg = pdg.GetParticle(ipdg);
 596     info->fEventNr = fEventNr;
 597     info->Update();
 598     //////////////////////////////////////////////////////////////////////
 599     br->SetAddress(&info);
 600     fTreeGenTracks->Fill();
 601   }
 602   fTreeGenTracks->AutoSave();
 603   if (fDebug > 2) cerr<<"end of TreeKLoop"<<endl;
 604   return 0;
 605 }
 606
 607
 608
 609
 610 ////////////////////////////////////////////////////////////////////////
 611 Int_t AliGenInfoMaker::TreeDLoop()
 612 {
 613   //
 614   // open the file with treeD
 615   // loop over all entries there and save information about some tracks
 616   //
 617
 618   Int_t nInnerSector = fParamTPC->GetNInnerSector();
 619   Int_t rowShift = 0;
 620   Int_t zero=fParamTPC->GetZeroSup()+6;
 621   //
 622   //
 623   AliSimDigits digitsAddress, *digits=&digitsAddress;
 624   fTreeD->GetBranch("Segment")->SetAddress(&digits);
 625
 626   Int_t sectorsByRows=(Int_t)fTreeD->GetEntries();
 627   if (fDebug > 1) cout<<"\tsectorsByRows = "<<sectorsByRows<<endl;
 628   for (Int_t i=0; i<sectorsByRows; i++) {
 629     if (!fTreeD->GetEvent(i)) continue;
 630     Int_t sec,row;
 631     fParamTPC->AdjustSectorRow(digits->GetID(),sec,row);
 632     if (fDebug > 1) cout<<sec<<' '<<row<<"                          \r";
 633     // here I expect that upper sectors follow lower sectors
 634     if (sec > nInnerSector) rowShift = fParamTPC->GetNRowLow();
 635     //
 636     digits->ExpandTrackBuffer();
 637     digits->First();
 638     do {
 639       Int_t iRow=digits->CurrentRow();
 640       Int_t iColumn=digits->CurrentColumn();
 641       Short_t digitValue = digits->CurrentDigit();
 642       if (digitValue >= zero) {
 643         Int_t label;
 644         for (Int_t j = 0; j<3; j++) {
 645           //      label = digits->GetTrackID(iRow,iColumn,j);
 646           label = digits->GetTrackIDFast(iRow,iColumn,j)-2;
 647           if (label >= (Int_t)fNParticles) { //don't label from bakground event
 648             continue;
 649           }
 650           if (label >= 0 && label <= (Int_t)fNParticles) {
 651             if (fDebug > 6 ) {
 652               cout<<"Inner loop: sector, iRow, iColumn, label, value, row "
 653                   <<sec<<" "
 654                   <<iRow<<" "<<iColumn<<" "<<label<<" "<<digitValue
 655                   <<" "<<row<<endl;
 656             }
 657             AliMCInfo * info = GetInfo(label);
 658             if (info){
 659               info->fTPCRow.SetRow(row+rowShift);
 660             }
 661           }
 662         }
 663       }
 664     } while (digits->Next());
 665   }
 666
 667   if (fDebug > 2) cerr<<"end of TreeDLoop"<<endl;
 668   return 0;
 669 }
 670
 671
 672 ////////////////////////////////////////////////////////////////////////
 673 Int_t AliGenInfoMaker::TreeTRLoop()
 674 {
 675   //
 676   // loop over TrackReferences and store the first one for each track
 677   //
 678   TTree * treeTR = fTreeTR;
 679   Int_t nPrimaries = (Int_t) treeTR->GetEntries();
 680   if (fDebug > 1) cout<<"There are "<<nPrimaries<<" entries in TreeTR"<<endl;
 681   //
 682   //
 683   //track references for TPC
 684   TClonesArray* TPCArrayTR = new TClonesArray("AliTrackReference");
 685   TClonesArray* ITSArrayTR = new TClonesArray("AliTrackReference");
 686   TClonesArray* TRDArrayTR = new TClonesArray("AliTrackReference");
 687   TClonesArray* RunArrayTR = new TClonesArray("AliTrackReference");
 688   //
 689   if (treeTR->GetBranch("TPC"))    treeTR->GetBranch("TPC")->SetAddress(&TPCArrayTR);
 690   if (treeTR->GetBranch("ITS"))    treeTR->GetBranch("ITS")->SetAddress(&ITSArrayTR);
 691   if (treeTR->GetBranch("TRD"))    treeTR->GetBranch("TRD")->SetAddress(&TRDArrayTR);
 692   if (treeTR->GetBranch("AliRun")) treeTR->GetBranch("AliRun")->SetAddress(&RunArrayTR);
 693   //
 694   //
 695   //
 696   for (Int_t iPrimPart = 0; iPrimPart<nPrimaries; iPrimPart++) {
 697     treeTR->GetEntry(iPrimPart);
 698     //
 699     // Loop over TPC references
 700     //
 701     for (Int_t iTrackRef = 0; iTrackRef < TPCArrayTR->GetEntriesFast(); iTrackRef++) {
 702       AliTrackReference *trackRef = (AliTrackReference*)TPCArrayTR->At(iTrackRef);
 703       //
 704       if (trackRef->Pt()<fgTPCPtCut) continue;
 705       Int_t label = trackRef->GetTrack();
 706       AliMCInfo * info = GetInfo(label);
 707       if (!info) info = MakeInfo(label);
 708       if (!info) continue;
 709       TClonesArray & arr = *(info->fTPCReferences);
 710       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 711     }
 712     //
 713     // Loop over ITS references
 714     //
 715     for (Int_t iTrackRef = 0; iTrackRef < ITSArrayTR->GetEntriesFast(); iTrackRef++) {
 716       AliTrackReference *trackRef = (AliTrackReference*)ITSArrayTR->At(iTrackRef);
 717       //
 718       if (trackRef->Pt()<fgTPCPtCut) continue;
 719       Int_t label = trackRef->GetTrack();
 720       AliMCInfo * info = GetInfo(label);
 721       if ( (!info) && trackRef->Pt()<fgITSPtCut) continue;
 722       if (!info) info = MakeInfo(label);
 723       if (!info) continue;
 724       TClonesArray & arr = *(info->fITSReferences);
 725       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 726     }
 727     //
 728     // Loop over TRD references
 729     //
 730     for (Int_t iTrackRef = 0; iTrackRef < TRDArrayTR->GetEntriesFast(); iTrackRef++) {
 731       AliTrackReference *trackRef = (AliTrackReference*)TRDArrayTR->At(iTrackRef);
 732       //
 733       if (trackRef->Pt()<fgTPCPtCut) continue;
 734       Int_t label = trackRef->GetTrack();
 735       AliMCInfo * info = GetInfo(label);
 736       if ( (!info) && trackRef->Pt()<fgTRDPtCut) continue;
 737       if (!info) info = MakeInfo(label);
 738       if (!info) continue;
 739       TClonesArray & arr = *(info->fTRDReferences);
 740       new (arr[arr.GetEntriesFast()]) AliTrackReference(*trackRef);
 741     }
 742     //
 743     // get dacay position
 744     for (Int_t iTrackRef = 0; iTrackRef < RunArrayTR->GetEntriesFast(); iTrackRef++) {
 745       AliTrackReference *trackRef = (AliTrackReference*)RunArrayTR->At(iTrackRef);
 746       //
 747       Int_t label = trackRef->GetTrack();
 748       AliMCInfo * info = GetInfo(label);
 749       if (!info) continue;
 750       if (!trackRef->TestBit(BIT(2))) continue;  //if not decay
 751       //      if (TMath::Abs(trackRef.X());
 752       info->fTRdecay = *trackRef;
 753     }
 754   }
 755   //
 756   TPCArrayTR->Delete();
 757   delete  TPCArrayTR;
 758   TRDArrayTR->Delete();
 759   delete  TRDArrayTR;
 760   ITSArrayTR->Delete();
 761   delete  ITSArrayTR;
 762   RunArrayTR->Delete();
 763   delete  RunArrayTR;
 764   //
 765   return 0;
 766 }
 767
 768 ////////////////////////////////////////////////////////////////////////
 769 Float_t AliGenInfoMaker::TR2LocalX(AliTrackReference *trackRef,
 770                                     AliTPCParam *paramTPC) {
 771
 772   Float_t x[3] = { trackRef->X(),trackRef->Y(),trackRef->Z()};
 773   Int_t index[4];
 774   paramTPC->Transform0to1(x,index);
 775   paramTPC->Transform1to2(x,index);
 776   return x[0];
 777 }
 778 ////////////////////////////////////////////////////////////////////////
 779
 780
 781
 782 TH1F * AliComparisonDraw::DrawXY(const char * chx, const char *chy, const char* selection,
 783                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
 784 {
 785   //
 786   Double_t* bins = CreateLogBins(nbins, minx, maxx);
 787   Int_t nBinsRes = 30;
 788   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, minx, maxx, nBinsRes, miny, maxy);
 789   char cut[1000];
 790   sprintf(cut,"%s&&%s",selection,quality);
 791   char expression[1000];
 792   sprintf(expression,"%s:%s>>hRes2",chy,chx);
 793   fTree->Draw(expression, cut, "groff");
 794   TH1F* hMean=0;
 795   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 796   AliLabelAxes(hRes, chx, chy);
 797   //
 798   delete hRes2;
 799   delete[] bins;
 800   fRes  = hRes;
 801   fMean = hMean;
 802   return hRes;
 803 }
 804
 805 TH1F * AliComparisonDraw::DrawLogXY(const char * chx, const char *chy, const char* selection,
 806                 const char * quality, Int_t nbins, Float_t minx, Float_t maxx, Float_t miny, Float_t maxy)
 807 {
 808   //
 809   Double_t* bins = CreateLogBins(nbins, minx, maxx);
 810   Int_t nBinsRes = 30;
 811   TH2F* hRes2 = new TH2F("hRes2", "residuals", nbins, bins, nBinsRes, miny, maxy);
 812   char cut[1000];
 813   sprintf(cut,"%s&&%s",selection,quality);
 814   char expression[1000];
 815   sprintf(expression,"%s:%s>>hRes2",chy,chx);
 816   fTree->Draw(expression, cut, "groff");
 817   TH1F* hMean=0;
 818   TH1F* hRes = CreateResHisto(hRes2, &hMean);
 819   AliLabelAxes(hRes, chx, chy);
 820   //
 821   delete hRes2;
 822   delete[] bins;
 823   fRes  = hRes;
 824   fMean = hMean;
 825   return hRes;
 826 }
 827
 828 ///////////////////////////////////////////////////////////////////////////////////
 829 ///////////////////////////////////////////////////////////////////////////////////
 830 TH1F * AliComparisonDraw::Eff(const char *variable, const char* selection, const char * quality,
 831                               Int_t nbins, Float_t min, Float_t max)
 832 {
 833   //
 834   //
 835   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, min, max);
 836   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, min, max);
 837   char inputGen[1000];
 838   sprintf(inputGen,"%s>>hGen", variable);
 839   fTree->Draw(inputGen, selection, "groff");
 840   char selectionRec[256];
 841   sprintf(selectionRec, "%s && %s", selection, quality);
 842   char inputRec[1000];
 843   sprintf(inputRec,"%s>>hRec", variable);
 844   fTree->Draw(inputRec, selectionRec, "groff");
 845   //
 846   TH1F* hEff = CreateEffHisto(hGen, hRec);
 847   AliLabelAxes(hEff, variable, "#epsilon [%]");
 848   fRes = hEff;
 849   delete hRec;
 850   delete hGen;
 851   return hEff;
 852 }
 853
 854
 855
 856 ///////////////////////////////////////////////////////////////////////////////////
 857 ///////////////////////////////////////////////////////////////////////////////////
 858 TH1F * AliComparisonDraw::EffLog(const char *variable, const char* selection, const char * quality,
 859                               Int_t nbins, Float_t min, Float_t max)
 860 {
 861   //
 862   //
 863   Double_t* bins = CreateLogBins(nbins, min, max);
 864   TH1F* hGen = new TH1F("hGen", "gen. tracks", nbins, bins);
 865   TH1F* hRec = new TH1F("hRec", "rec. tracks", nbins, bins);
 866   char inputGen[1000];
 867   sprintf(inputGen,"%s>>hGen", variable);
 868   fTree->Draw(inputGen, selection, "groff");
 869   char selectionRec[256];
 870   sprintf(selectionRec, "%s && %s", selection, quality);
 871   char inputRec[1000];
 872   sprintf(inputRec,"%s>>hRec", variable);
 873   fTree->Draw(inputRec, selectionRec, "groff");
 874   //
 875   TH1F* hEff = CreateEffHisto(hGen, hRec);
 876   AliLabelAxes(hEff, variable, "#epsilon [%]");
 877   fRes = hEff;
 878   delete hRec;
 879   delete hGen;
 880   delete[] bins;
 881   return hEff;
 882 }
 883
 884
 885 ///////////////////////////////////////////////////////////////////////////////////
 886 ///////////////////////////////////////////////////////////////////////////////////
 887
 888 Double_t* AliComparisonDraw::CreateLogBins(Int_t nBins, Double_t xMin, Double_t xMax)
 889 {
 890   Double_t* bins = new Double_t[nBins+1];
 891   bins[0] = xMin;
 892   Double_t factor = pow(xMax/xMin, 1./nBins);
 893   for (Int_t i = 1; i <= nBins; i++)
 894     bins[i] = factor * bins[i-1];
 895   return bins;
 896 }
 897
 898
 899
 900
 901 void AliComparisonDraw::AliLabelAxes(TH1* histo, const char* xAxisTitle, const char* yAxisTitle)
 902 {
 903   //
 904   histo->GetXaxis()->SetTitle(xAxisTitle);
 905   histo->GetYaxis()->SetTitle(yAxisTitle);
 906 }
 907
 908
 909 TH1F* AliComparisonDraw::CreateEffHisto(TH1F* hGen, TH1F* hRec)
 910 {
 911   //
 912   Int_t nBins = hGen->GetNbinsX();
 913   TH1F* hEff = (TH1F*) hGen->Clone("hEff");
 914   hEff->SetTitle("");
 915   hEff->SetStats(kFALSE);
 916   hEff->SetMinimum(0.);
 917   hEff->SetMaximum(110.);
 918   //
 919   for (Int_t iBin = 0; iBin <= nBins; iBin++) {
 920     Double_t nGen = hGen->GetBinContent(iBin);
 921     Double_t nRec = hRec->GetBinContent(iBin);
 922     if (nGen > 0) {
 923       Double_t eff = nRec/nGen;
 924       hEff->SetBinContent(iBin, 100. * eff);
 925       Double_t error = sqrt((eff*(1.-eff)+0.01) / nGen);
 926       //      if (error == 0) error = sqrt(0.1/nGen);
 927       //
 928       if (error == 0) error = 0.0001;
 929       hEff->SetBinError(iBin, 100. * error);
 930     } else {
 931       hEff->SetBinContent(iBin, 100. * 0.5);
 932       hEff->SetBinError(iBin, 100. * 0.5);
 933     }
 934   }
 935   return hEff;
 936 }
 937
 938
 939
 940 TH1F* AliComparisonDraw::CreateResHisto(TH2F* hRes2, TH1F **phMean,  Bool_t drawBinFits,
 941                      Bool_t overflowBinFits)
 942 {
 943   TVirtualPad* currentPad = gPad;
 944   TAxis* axis = hRes2->GetXaxis();
 945   Int_t nBins = axis->GetNbins();
 946   TH1F* hRes, *hMean;
 947   if (axis->GetXbins()->GetSize()){
 948     hRes = new TH1F("hRes", "", nBins, axis->GetXbins()->GetArray());
 949     hMean = new TH1F("hMean", "", nBins, axis->GetXbins()->GetArray());
 950   }
 951   else{
 952     hRes = new TH1F("hRes", "", nBins, axis->GetXmin(), axis->GetXmax());
 953     hMean = new TH1F("hMean", "", nBins, axis->GetXmin(), axis->GetXmax());
 954
 955   }
 956   hRes->SetStats(false);
 957   hRes->SetOption("E");
 958   hRes->SetMinimum(0.);
 959   //
 960   hMean->SetStats(false);
 961   hMean->SetOption("E");
 962
 963   // create the fit function
 964   //TKFitGaus* fitFunc = new TKFitGaus("resFunc");
 965   //   TF1 * fitFunc = new TF1("G","[3]+[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 966   TF1 * fitFunc = new TF1("G","[0]*exp(-(x-[1])*(x-[1])/(2.*[2]*[2]))",-3,3);
 967
 968   fitFunc->SetLineWidth(2);
 969   fitFunc->SetFillStyle(0);
 970   // create canvas for fits
 971   TCanvas* canBinFits = NULL;
 972   Int_t nPads = (overflowBinFits) ? nBins+2 : nBins;
 973   Int_t nx = Int_t(sqrt(nPads-1.));// + 1;
 974   Int_t ny = (nPads-1) / nx + 1;
 975   if (drawBinFits) {
 976     canBinFits = (TCanvas*)gROOT->FindObject("canBinFits");
 977     if (canBinFits) delete canBinFits;
 978     canBinFits = new TCanvas("canBinFits", "fits of bins", 200, 100, 500, 700);
 979     canBinFits->Divide(nx, ny);
 980   }
 981
 982   // loop over x bins and fit projection
 983   Int_t dBin = ((overflowBinFits) ? 1 : 0);
 984   for (Int_t bin = 1-dBin; bin <= nBins+dBin; bin++) {
 985     if (drawBinFits) canBinFits->cd(bin + dBin);
 986     TH1D* hBin = hRes2->ProjectionY("hBin", bin, bin);
 987     //
 988     if (hBin->GetEntries() > 5) {
 989       fitFunc->SetParameters(hBin->GetMaximum(),hBin->GetMean(),hBin->GetRMS());
 990       hBin->Fit(fitFunc,"s");
 991       Double_t sigma = TMath::Abs(fitFunc->GetParameter(2));
 992
 993       if (sigma > 0.){
 994         hRes->SetBinContent(bin, TMath::Abs(fitFunc->GetParameter(2)));
 995         hMean->SetBinContent(bin, fitFunc->GetParameter(1));
 996       }
 997       else{
 998         hRes->SetBinContent(bin, 0.);
 999         hMean->SetBinContent(bin,0);
1000       }
1001       hRes->SetBinError(bin, fitFunc->GetParError(2));
1002       hMean->SetBinError(bin, fitFunc->GetParError(1));
1003
1004       //
1005       //
1006
1007     } else {
1008       hRes->SetBinContent(bin, 0.);
1009       hRes->SetBinError(bin, 0.);
1010       hMean->SetBinContent(bin, 0.);
1011       hMean->SetBinError(bin, 0.);
1012     }
1013
1014
1015     if (drawBinFits) {
1016       char name[256];
1017       if (bin == 0) {
1018         sprintf(name, "%s < %.4g", axis->GetTitle(), axis->GetBinUpEdge(bin));
1019       } else if (bin == nBins+1) {
1020         sprintf(name, "%.4g < %s", axis->GetBinLowEdge(bin), axis->GetTitle());
1021       } else {
1022         sprintf(name, "%.4g < %s < %.4g", axis->GetBinLowEdge(bin),
1023                 axis->GetTitle(), axis->GetBinUpEdge(bin));
1024       }
1025       canBinFits->cd(bin + dBin);
1026       hBin->SetTitle(name);
1027       hBin->SetStats(kTRUE);
1028       hBin->DrawCopy("E");
1029       canBinFits->Update();
1030       canBinFits->Modified();
1031       canBinFits->Update();
1032     }
1033
1034     delete hBin;
1035   }
1036
1037   delete fitFunc;
1038   currentPad->cd();
1039   *phMean = hMean;
1040   return hRes;
1041 }
1042
1043
1044