TPC/AliTPC.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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Time Projection Chamber                                                  //
  21 //  This class contains the basic functions for the Time Projection Chamber  //
  22 //  detector. Functions specific to one particular geometry are              //
  23 //  contained in the derived classes                                         //
  24 //                                                                           //
  25 //Begin_Html
  26 /*
  27 <img src="picts/AliTPCClass.gif">
  28 */
  29 //End_Html
  30 //                                                                           //
  31 //                                                                          //
  32 ///////////////////////////////////////////////////////////////////////////////
  33
  34 //
  35
  36 #include <Riostream.h>
  37 #include <stdlib.h>
  38
  39 #include <TFile.h>
  40 #include <TGeometry.h>
  41 #include <TInterpreter.h>
  42 #include <TMath.h>
  43 #include <TMatrix.h>
  44 #include <TNode.h>
  45 #include <TObjectTable.h>
  46 #include <TParticle.h>
  47 #include <TROOT.h>
  48 #include <TRandom.h>
  49 #include <TSystem.h>
  50 #include <TTUBS.h>
  51 #include <TTree.h>
  52 #include <TVector.h>
  53 #include <TVirtualMC.h>
  54 #include <TString.h>
  55 #include <TF2.h>
  56
  57 #include "AliArrayBranch.h"
  58 #include "AliClusters.h"
  59 #include "AliComplexCluster.h"
  60 #include "AliDigits.h"
  61 #include "AliMagF.h"
  62 #include "AliPoints.h"
  63 #include "AliRun.h"
  64 #include "AliRunLoader.h"
  65 #include "AliSimDigits.h"
  66 #include "AliTPC.h"
  67 #include "AliTPC.h"
  68 #include "AliTPCClustersArray.h"
  69 #include "AliTPCClustersRow.h"
  70 #include "AliTPCDigitsArray.h"
  71 #include "AliTPCLoader.h"
  72 #include "AliTPCPRF2D.h"
  73 #include "AliTPCParamSR.h"
  74 #include "AliTPCRF1D.h"
  75 #include "AliTPCTrackHits.h"
  76 #include "AliTPCTrackHitsV2.h"
  77 #include "AliTPCcluster.h"
  78 #include "AliTrackReference.h"
  79 #include "AliMC.h"
  80 #include "AliTPCDigitizer.h"
  81
  82
  83 ClassImp(AliTPC)
  84
  85 //_____________________________________________________________________________
  86 // helper class for fast matrix and vector manipulation - no range checking
  87 // origin - Marian Ivanov
  88
  89 class AliTPCFastMatrix : public TMatrix {
  90 public :
  91   AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb);
  92   Float_t & UncheckedAt(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces
  93   Float_t   UncheckedAtFast(Int_t rown, Int_t coln) const  {return  (fIndex[coln])[rown];} //fast acces
  94 };
  95
  96 AliTPCFastMatrix::AliTPCFastMatrix(Int_t rowlwb, Int_t rowupb, Int_t collwb, Int_t colupb):
  97   TMatrix(rowlwb, rowupb,collwb,colupb)
  98    {
  99    };
 100 //_____________________________________________________________________________
 101 class AliTPCFastVector : public TVector {
 102 public :
 103   AliTPCFastVector(Int_t size);
 104   virtual ~AliTPCFastVector(){;}
 105   Float_t & UncheckedAt(Int_t index) const  {return  fElements[index];} //fast acces
 106
 107 };
 108
 109 AliTPCFastVector::AliTPCFastVector(Int_t size):TVector(size){
 110 };
 111
 112 //_____________________________________________________________________________
 113 AliTPC::AliTPC()
 114 {
 115   //
 116   // Default constructor
 117   //
 118   fIshunt   = 0;
 119   fHits     = 0;
 120   fDigits   = 0;
 121   fNsectors = 0;
 122   fDigitsArray = 0;
 123   fClustersArray = 0;
 124   fDefaults = 0;
 125   fTrackHits = 0;
 126   fTrackHitsOld = 0;
 127   fHitType = 2; //default CONTAINERS - based on ROOT structure
 128   fTPCParam = 0;
 129   fNoiseTable = 0;
 130   fActiveSectors =0;
 131
 132 }
 133
 134 //_____________________________________________________________________________
 135 AliTPC::AliTPC(const char *name, const char *title)
 136       : AliDetector(name,title)
 137 {
 138   //
 139   // Standard constructor
 140   //
 141
 142   //
 143   // Initialise arrays of hits and digits
 144   fHits     = new TClonesArray("AliTPChit",  176);
 145   gAlice->GetMCApp()->AddHitList(fHits);
 146   fDigitsArray = 0;
 147   fClustersArray= 0;
 148   fDefaults = 0;
 149   //
 150   fTrackHits = new AliTPCTrackHitsV2;
 151   fTrackHits->SetHitPrecision(0.002);
 152   fTrackHits->SetStepPrecision(0.003);
 153   fTrackHits->SetMaxDistance(100);
 154
 155   fTrackHitsOld = new AliTPCTrackHits;  //MI - 13.09.2000
 156   fTrackHitsOld->SetHitPrecision(0.002);
 157   fTrackHitsOld->SetStepPrecision(0.003);
 158   fTrackHitsOld->SetMaxDistance(100);
 159
 160   fNoiseTable =0;
 161
 162   fHitType = 2;
 163   fActiveSectors = 0;
 164   //
 165   // Initialise counters
 166   fNsectors = 0;
 167
 168   //
 169   fIshunt     =  0;
 170   //
 171   // Initialise color attributes
 172   SetMarkerColor(kYellow);
 173
 174   //
 175   //  Set TPC parameters
 176   //
 177
 178
 179   if (!strcmp(title,"Default")) {
 180     fTPCParam = new AliTPCParamSR;
 181   } else {
 182     cerr<<"AliTPC warning: in Config.C you must set non-default parameters\n";
 183     fTPCParam=0;
 184   }
 185
 186 }
 187
 188 //_____________________________________________________________________________
 189 AliTPC::AliTPC(const AliTPC& t):AliDetector(t){
 190   //
 191   // dummy copy constructor
 192   //
 193 }
 194 AliTPC::~AliTPC()
 195 {
 196   //
 197   // TPC destructor
 198   //
 199
 200   fIshunt   = 0;
 201   delete fHits;
 202   delete fDigits;
 203   delete fTPCParam;
 204   delete fTrackHits; //MI 15.09.2000
 205   delete fTrackHitsOld; //MI 10.12.2001
 206   if (fNoiseTable) delete [] fNoiseTable;
 207
 208 }
 209
 210 //_____________________________________________________________________________
 211 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 212 {
 213   //
 214   // Add a hit to the list
 215   //
 216   //  TClonesArray &lhits = *fHits;
 217   //  new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 218   if (fHitType&1){
 219     TClonesArray &lhits = *fHits;
 220     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 221   }
 222   if (fHitType>1)
 223    AddHit2(track,vol,hits);
 224 }
 225
 226 //_____________________________________________________________________________
 227 void AliTPC::BuildGeometry()
 228 {
 229
 230   //
 231   // Build TPC ROOT TNode geometry for the event display
 232   //
 233   TNode *nNode, *nTop;
 234   TTUBS *tubs;
 235   Int_t i;
 236   const int kColorTPC=19;
 237   char name[5], title[25];
 238   const Double_t kDegrad=TMath::Pi()/180;
 239   const Double_t kRaddeg=180./TMath::Pi();
 240
 241
 242   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
 243   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
 244
 245   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
 246   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
 247
 248   Int_t nLo = fTPCParam->GetNInnerSector()/2;
 249   Int_t nHi = fTPCParam->GetNOuterSector()/2;
 250
 251   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
 252   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
 253   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
 254   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
 255
 256
 257   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
 258   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
 259
 260   Double_t rl,ru;
 261
 262
 263   //
 264   // Get ALICE top node
 265   //
 266
 267   nTop=gAlice->GetGeometry()->GetNode("alice");
 268
 269   //  inner sectors
 270
 271   rl = fTPCParam->GetInnerRadiusLow();
 272   ru = fTPCParam->GetInnerRadiusUp();
 273
 274
 275   for(i=0;i<nLo;i++) {
 276     sprintf(name,"LS%2.2d",i);
 277     name[4]='\0';
 278     sprintf(title,"TPC low sector %3d",i);
 279     title[24]='\0';
 280
 281     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
 282                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
 283     tubs->SetNumberOfDivisions(1);
 284     nTop->cd();
 285     nNode = new TNode(name,title,name,0,0,0,"");
 286     nNode->SetLineColor(kColorTPC);
 287     fNodes->Add(nNode);
 288   }
 289
 290   // Outer sectors
 291
 292   rl = fTPCParam->GetOuterRadiusLow();
 293   ru = fTPCParam->GetOuterRadiusUp();
 294
 295   for(i=0;i<nHi;i++) {
 296     sprintf(name,"US%2.2d",i);
 297     name[4]='\0';
 298     sprintf(title,"TPC upper sector %d",i);
 299     title[24]='\0';
 300     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
 301                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
 302     tubs->SetNumberOfDivisions(1);
 303     nTop->cd();
 304     nNode = new TNode(name,title,name,0,0,0,"");
 305     nNode->SetLineColor(kColorTPC);
 306     fNodes->Add(nNode);
 307   }
 308
 309 }
 310
 311 //_____________________________________________________________________________
 312 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) const
 313 {
 314   //
 315   // Calculate distance from TPC to mouse on the display
 316   // Dummy procedure
 317   //
 318   return 9999;
 319 }
 320
 321 void AliTPC::Clusters2Tracks() const
 322  {
 323   //-----------------------------------------------------------------
 324   // This is a track finder.
 325   //-----------------------------------------------------------------
 326   Error("Clusters2Tracks",
 327   "Dummy function !  Call AliTPCtracker::Clusters2Tracks(...) instead !");
 328  }
 329
 330 //_____________________________________________________________________________
 331 void AliTPC::CreateMaterials()
 332 {
 333   //-----------------------------------------------
 334   // Create Materials for for TPC simulations
 335   //-----------------------------------------------
 336
 337   //-----------------------------------------------------------------
 338   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 339   //-----------------------------------------------------------------
 340
 341   Int_t iSXFLD=gAlice->Field()->Integ();
 342   Float_t sXMGMX=gAlice->Field()->Max();
 343
 344   Float_t amat[5]; // atomic numbers
 345   Float_t zmat[5]; // z
 346   Float_t wmat[5]; // proportions
 347
 348   Float_t density;
 349   Float_t apure[2];
 350
 351
 352   //***************** Gases *************************
 353
 354   //-------------------------------------------------
 355   // pure gases
 356   //-------------------------------------------------
 357
 358   // Neon
 359
 360
 361   amat[0]= 20.18;
 362   zmat[0]= 10.;
 363   density = 0.0009;
 364
 365   apure[0]=amat[0];
 366
 367   AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
 368
 369   // Argon
 370
 371   amat[0]= 39.948;
 372   zmat[0]= 18.;
 373   density = 0.001782;
 374
 375   apure[1]=amat[0];
 376
 377   AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
 378
 379
 380   //--------------------------------------------------------------
 381   // gases - compounds
 382   //--------------------------------------------------------------
 383
 384   Float_t amol[3];
 385
 386   // CO2
 387
 388   amat[0]=12.011;
 389   amat[1]=15.9994;
 390
 391   zmat[0]=6.;
 392   zmat[1]=8.;
 393
 394   wmat[0]=1.;
 395   wmat[1]=2.;
 396
 397   density=0.001977;
 398
 399   amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
 400
 401   AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
 402
 403   // CF4
 404
 405   amat[0]=12.011;
 406   amat[1]=18.998;
 407
 408   zmat[0]=6.;
 409   zmat[1]=9.;
 410
 411   wmat[0]=1.;
 412   wmat[1]=4.;
 413
 414   density=0.003034;
 415
 416   amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
 417
 418   AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
 419
 420
 421   // CH4
 422
 423   amat[0]=12.011;
 424   amat[1]=1.;
 425
 426   zmat[0]=6.;
 427   zmat[1]=1.;
 428
 429   wmat[0]=1.;
 430   wmat[1]=4.;
 431
 432   density=0.000717;
 433
 434   amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
 435
 436   AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
 437
 438   //----------------------------------------------------------------
 439   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
 440   //----------------------------------------------------------------
 441
 442   char namate[21]="";
 443   density = 0.;
 444   Float_t am=0;
 445   Int_t nc;
 446   Float_t rho,absl,x0,buf[1];
 447   Int_t nbuf;
 448   Float_t a,z;
 449
 450   for(nc = 0;nc<fNoComp;nc++)
 451     {
 452
 453       // retrive material constants
 454
 455       gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,x0,absl,buf,nbuf);
 456
 457       amat[nc] = a;
 458       zmat[nc] = z;
 459
 460       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 461
 462       am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
 463       density += fMixtProp[nc]*rho;  // density of the mixture
 464
 465     }
 466
 467   // mixture proportions by weight!
 468
 469   for(nc = 0;nc<fNoComp;nc++)
 470     {
 471
 472       Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 473
 474       wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
 475                  apure[nnc] : amol[nnc])/am;
 476
 477     }
 478
 479   // Drift gases 1 - nonsensitive, 2 - sensitive
 480
 481   AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
 482   AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
 483
 484
 485   // Air
 486
 487   amat[0] = 14.61;
 488   zmat[0] = 7.3;
 489   density = 0.001205;
 490
 491   AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
 492
 493
 494   //----------------------------------------------------------------------
 495   //               solid materials
 496   //----------------------------------------------------------------------
 497
 498
 499   // Kevlar C14H22O2N2
 500
 501   amat[0] = 12.011;
 502   amat[1] = 1.;
 503   amat[2] = 15.999;
 504   amat[3] = 14.006;
 505
 506   zmat[0] = 6.;
 507   zmat[1] = 1.;
 508   zmat[2] = 8.;
 509   zmat[3] = 7.;
 510
 511   wmat[0] = 14.;
 512   wmat[1] = 22.;
 513   wmat[2] = 2.;
 514   wmat[3] = 2.;
 515
 516   density = 1.45;
 517
 518   AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
 519
 520   // NOMEX
 521
 522   amat[0] = 12.011;
 523   amat[1] = 1.;
 524   amat[2] = 15.999;
 525   amat[3] = 14.006;
 526
 527   zmat[0] = 6.;
 528   zmat[1] = 1.;
 529   zmat[2] = 8.;
 530   zmat[3] = 7.;
 531
 532   wmat[0] = 14.;
 533   wmat[1] = 22.;
 534   wmat[2] = 2.;
 535   wmat[3] = 2.;
 536
 537   density = 0.03;
 538
 539
 540   AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
 541
 542   // Makrolon C16H18O3
 543
 544   amat[0] = 12.011;
 545   amat[1] = 1.;
 546   amat[2] = 15.999;
 547
 548   zmat[0] = 6.;
 549   zmat[1] = 1.;
 550   zmat[2] = 8.;
 551
 552   wmat[0] = 16.;
 553   wmat[1] = 18.;
 554   wmat[2] = 3.;
 555
 556   density = 1.2;
 557
 558   AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
 559
 560   // Mylar C5H4O2
 561
 562   amat[0]=12.011;
 563   amat[1]=1.;
 564   amat[2]=15.9994;
 565
 566   zmat[0]=6.;
 567   zmat[1]=1.;
 568   zmat[2]=8.;
 569
 570   wmat[0]=5.;
 571   wmat[1]=4.;
 572   wmat[2]=2.;
 573
 574   density = 1.39;
 575
 576   AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
 577
 578   // SiO2 - used later for the glass fiber
 579
 580   amat[0]=28.086;
 581   amat[1]=15.9994;
 582
 583   zmat[0]=14.;
 584   zmat[1]=8.;
 585
 586   wmat[0]=1.;
 587   wmat[1]=2.;
 588
 589
 590   AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2)
 591
 592   // Al
 593
 594   amat[0] = 26.98;
 595   zmat[0] = 13.;
 596
 597   density = 2.7;
 598
 599   AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
 600
 601   // Si
 602
 603   amat[0] = 28.086;
 604   zmat[0] = 14.;
 605
 606   density = 2.33;
 607
 608   AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
 609
 610   // Cu
 611
 612   amat[0] = 63.546;
 613   zmat[0] = 29.;
 614
 615   density = 8.96;
 616
 617   AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
 618
 619   // Tedlar C2H3F
 620
 621   amat[0] = 12.011;
 622   amat[1] = 1.;
 623   amat[2] = 18.998;
 624
 625   zmat[0] = 6.;
 626   zmat[1] = 1.;
 627   zmat[2] = 9.;
 628
 629   wmat[0] = 2.;
 630   wmat[1] = 3.;
 631   wmat[2] = 1.;
 632
 633   density = 1.71;
 634
 635   AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
 636
 637
 638   // Plexiglas  C5H8O2
 639
 640   amat[0]=12.011;
 641   amat[1]=1.;
 642   amat[2]=15.9994;
 643
 644   zmat[0]=6.;
 645   zmat[1]=1.;
 646   zmat[2]=8.;
 647
 648   wmat[0]=5.;
 649   wmat[1]=8.;
 650   wmat[2]=2.;
 651
 652   density=1.18;
 653
 654   AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
 655
 656   // Epoxy - C14 H20 O3
 657
 658
 659   amat[0]=12.011;
 660   amat[1]=1.;
 661   amat[2]=15.9994;
 662
 663   zmat[0]=6.;
 664   zmat[1]=1.;
 665   zmat[2]=8.;
 666
 667   wmat[0]=14.;
 668   wmat[1]=20.;
 669   wmat[2]=3.;
 670
 671   density=1.25;
 672
 673   AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat);
 674
 675   // Carbon
 676
 677   amat[0]=12.011;
 678   zmat[0]=6.;
 679   density= 2.265;
 680
 681   AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.);
 682
 683   // get epoxy
 684
 685   gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,x0,absl,buf,nbuf);
 686
 687   // Carbon fiber
 688
 689   wmat[0]=0.644; // by weight!
 690   wmat[1]=0.356;
 691
 692   density=0.5*(1.25+2.265);
 693
 694   AliMixture(47,"Cfiber",amat,zmat,density,2,wmat);
 695
 696   // get SiO2
 697
 698   gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,x0,absl,buf,nbuf);
 699
 700   wmat[0]=0.725; // by weight!
 701   wmat[1]=0.275;
 702
 703   density=1.7;
 704
 705   AliMixture(39,"G10",amat,zmat,density,2,wmat);
 706
 707
 708
 709
 710   //----------------------------------------------------------
 711   // tracking media for gases
 712   //----------------------------------------------------------
 713
 714   AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
 715   AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 716   AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 717   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
 718
 719   //-----------------------------------------------------------
 720   // tracking media for solids
 721   //-----------------------------------------------------------
 722
 723   AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 724   AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 725   AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 726   AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 727   AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 728   AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 729   AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 730   AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 731   AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 732   AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 733   AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 734   AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 735
 736 }
 737
 738 void AliTPC::GenerNoise(Int_t tablesize)
 739 {
 740   //
 741   //Generate table with noise
 742   //
 743   if (fTPCParam==0) {
 744     // error message
 745     fNoiseDepth=0;
 746     return;
 747   }
 748   if (fNoiseTable)  delete[] fNoiseTable;
 749   fNoiseTable = new Float_t[tablesize];
 750   fNoiseDepth = tablesize;
 751   fCurrentNoise =0; //!index of the noise in  the noise table
 752
 753   Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
 754   for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);
 755 }
 756
 757 Float_t AliTPC::GetNoise()
 758 {
 759   // get noise from table
 760   //  if ((fCurrentNoise%10)==0)
 761   //  fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
 762   if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
 763   return fNoiseTable[fCurrentNoise++];
 764   //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
 765 }
 766
 767
 768 Bool_t  AliTPC::IsSectorActive(Int_t sec) const
 769 {
 770   //
 771   // check if the sector is active
 772   if (!fActiveSectors) return kTRUE;
 773   else return fActiveSectors[sec];
 774 }
 775
 776 void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
 777 {
 778   // activate interesting sectors
 779   SetTreeAddress();//just for security
 780   if (fActiveSectors) delete [] fActiveSectors;
 781   fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 782   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 783   for (Int_t i=0;i<n;i++)
 784     if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector())  fActiveSectors[sectors[i]]=kTRUE;
 785
 786 }
 787
 788 void    AliTPC::SetActiveSectors(Int_t flag)
 789 {
 790   //
 791   // activate sectors which were hitted by tracks
 792   //loop over tracks
 793   SetTreeAddress();//just for security
 794   if (fHitType==0) return;  // if Clones hit - not short volume ID information
 795   if (fActiveSectors) delete [] fActiveSectors;
 796   fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 797   if (flag) {
 798     for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
 799     return;
 800   }
 801   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 802   TBranch * branch=0;
 803   if (TreeH() == 0x0)
 804    {
 805      Fatal("SetActiveSectors","Can not find TreeH in folder");
 806      return;
 807    }
 808   if (fHitType>1) branch = TreeH()->GetBranch("TPC2");
 809   else branch = TreeH()->GetBranch("TPC");
 810   Stat_t ntracks = TreeH()->GetEntries();
 811   // loop over all hits
 812   cout<<"\nAliTPC::SetActiveSectors():  Got "<<ntracks<<" tracks\n";
 813
 814   for(Int_t track=0;track<ntracks;track++)
 815    {
 816     ResetHits();
 817     //
 818     if (fTrackHits && fHitType&4) {
 819       TBranch * br1 = TreeH()->GetBranch("fVolumes");
 820       TBranch * br2 = TreeH()->GetBranch("fNVolumes");
 821       br1->GetEvent(track);
 822       br2->GetEvent(track);
 823       Int_t *volumes = fTrackHits->GetVolumes();
 824       for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++)
 825         fActiveSectors[volumes[j]]=kTRUE;
 826     }
 827
 828     //
 829     if (fTrackHitsOld && fHitType&2) {
 830       TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
 831       br->GetEvent(track);
 832       AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
 833       for (UInt_t j=0;j<ar->GetSize();j++){
 834         fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
 835       }
 836     }
 837   }
 838
 839 }
 840
 841
 842
 843
 844 void AliTPC::Digits2Clusters(Int_t /*eventnumber*/) const
 845 {
 846   //-----------------------------------------------------------------
 847   // This is a simple cluster finder.
 848   //-----------------------------------------------------------------
 849   Error("Digits2Clusters",
 850   "Dummy function !  Call AliTPCclusterer::Digits2Clusters(...) instead !");
 851 }
 852
 853 extern Double_t SigmaY2(Double_t, Double_t, Double_t);
 854 extern Double_t SigmaZ2(Double_t, Double_t);
 855 //_____________________________________________________________________________
 856 void AliTPC::Hits2Clusters(Int_t /*eventn*/)
 857 {
 858   //--------------------------------------------------------
 859   // TPC simple cluster generator from hits
 860   // obtained from the TPC Fast Simulator
 861   // The point errors are taken from the parametrization
 862   //--------------------------------------------------------
 863
 864   //-----------------------------------------------------------------
 865   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 866   //-----------------------------------------------------------------
 867   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
 868   // Jouri.Belikov@cern.ch
 869   //----------------------------------------------------------------
 870
 871   /////////////////////////////////////////////////////////////////////////////
 872   //
 873   //---------------------------------------------------------------------
 874   //   ALICE TPC Cluster Parameters
 875   //--------------------------------------------------------------------
 876
 877
 878
 879   // Cluster width in rphi
 880   const Float_t kACrphi=0.18322;
 881   const Float_t kBCrphi=0.59551e-3;
 882   const Float_t kCCrphi=0.60952e-1;
 883   // Cluster width in z
 884   const Float_t kACz=0.19081;
 885   const Float_t kBCz=0.55938e-3;
 886   const Float_t kCCz=0.30428;
 887
 888
 889   if (!fLoader) {
 890      cerr<<"AliTPC::Hits2Clusters(): output file not open !\n";
 891      return;
 892   }
 893
 894   //if(fDefaults == 0) SetDefaults();
 895
 896   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
 897   //
 898   TParticle *particle; // pointer to a given particle
 899   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 900   Int_t sector;
 901   Int_t ipart;
 902   Float_t xyz[5];
 903   Float_t pl,pt,tanth,rpad,ratio;
 904   Float_t cph,sph;
 905
 906   //---------------------------------------------------------------
 907   //  Get the access to the tracks
 908   //---------------------------------------------------------------
 909
 910   TTree *tH = TreeH();
 911   if (tH == 0x0)
 912    {
 913      Fatal("Hits2Clusters","Can not find TreeH in folder");
 914      return;
 915    }
 916   SetTreeAddress();
 917
 918   Stat_t ntracks = tH->GetEntries();
 919
 920   //Switch to the output file
 921
 922   if (fLoader->TreeR() == 0x0) fLoader->MakeTree("R");
 923
 924   cout<<"fTPCParam->GetTitle() = "<<fTPCParam->GetTitle()<<endl;
 925
 926   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
 927   rl->CdGAFile();
 928   //fTPCParam->Write(fTPCParam->GetTitle());
 929
 930   AliTPCClustersArray carray;
 931   carray.Setup(fTPCParam);
 932   carray.SetClusterType("AliTPCcluster");
 933   carray.MakeTree(fLoader->TreeR());
 934
 935   Int_t nclusters=0; //cluster counter
 936
 937   //------------------------------------------------------------
 938   // Loop over all sectors (72 sectors for 20 deg
 939   // segmentation for both lower and upper sectors)
 940   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
 941   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
 942   //
 943   // First cluster for sector 0 starts at "0"
 944   //------------------------------------------------------------
 945
 946   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++){
 947     //MI change
 948     fTPCParam->AdjustCosSin(isec,cph,sph);
 949
 950     //------------------------------------------------------------
 951     // Loop over tracks
 952     //------------------------------------------------------------
 953
 954     for(Int_t track=0;track<ntracks;track++){
 955       ResetHits();
 956       tH->GetEvent(track);
 957       //
 958       //  Get number of the TPC hits
 959       //
 960        tpcHit = (AliTPChit*)FirstHit(-1);
 961
 962       // Loop over hits
 963       //
 964        while(tpcHit){
 965
 966          if (tpcHit->fQ == 0.) {
 967            tpcHit = (AliTPChit*) NextHit();
 968            continue; //information about track (I.Belikov)
 969          }
 970         sector=tpcHit->fSector; // sector number
 971
 972        if(sector != isec){
 973          tpcHit = (AliTPChit*) NextHit();
 974          continue;
 975        }
 976         ipart=tpcHit->Track();
 977         particle=gAlice->GetMCApp()->Particle(ipart);
 978         pl=particle->Pz();
 979         pt=particle->Pt();
 980         if(pt < 1.e-9) pt=1.e-9;
 981         tanth=pl/pt;
 982         tanth = TMath::Abs(tanth);
 983         rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
 984         ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 985
 986         //   space-point resolutions
 987
 988         sigmaRphi=SigmaY2(rpad,tanth,pt);
 989         sigmaZ   =SigmaZ2(rpad,tanth   );
 990
 991         //   cluster widths
 992
 993         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
 994         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
 995
 996         // temporary protection
 997
 998         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
 999         if(sigmaZ < 0.) sigmaZ=0.4e-3;
1000         if(clRphi < 0.) clRphi=2.5e-3;
1001         if(clZ < 0.) clZ=2.5e-5;
1002
1003         //
1004
1005         //
1006         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
1007         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
1008         //
1009         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
1010         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
1011         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
1012           Float_t alpha=(isec < fTPCParam->GetNInnerSector()) ?
1013           fTPCParam->GetInnerAngle() : fTPCParam->GetOuterAngle();
1014           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
1015           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
1016         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
1017           if (TMath::Abs(xyz[1])>fTPCParam->GetZLength()) xyz[1]=tpcHit->Z();
1018         xyz[2]=sigmaRphi;                                     // fSigmaY2
1019         xyz[3]=sigmaZ;                                        // fSigmaZ2
1020         xyz[4]=tpcHit->fQ;                                    // q
1021
1022         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
1023         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
1024
1025         Int_t tracks[3]={tpcHit->Track(), -1, -1};
1026         AliTPCcluster cluster(tracks,xyz);
1027
1028         clrow->InsertCluster(&cluster); nclusters++;
1029
1030         tpcHit = (AliTPChit*)NextHit();
1031
1032
1033       } // end of loop over hits
1034
1035     }   // end of loop over tracks
1036
1037     Int_t nrows=fTPCParam->GetNRow(isec);
1038     for (Int_t irow=0; irow<nrows; irow++) {
1039         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
1040         if (!clrow) continue;
1041         carray.StoreRow(isec,irow);
1042         carray.ClearRow(isec,irow);
1043     }
1044
1045   } // end of loop over sectors
1046
1047   cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
1048   fLoader->WriteRecPoints("OVERWRITE");
1049
1050
1051 } // end of function
1052
1053 //_________________________________________________________________
1054 void AliTPC::Hits2ExactClustersSector(Int_t isec)
1055 {
1056   //--------------------------------------------------------
1057   //calculate exact cross point of track and given pad row
1058   //resulting values are expressed in "digit" coordinata
1059   //--------------------------------------------------------
1060
1061   //-----------------------------------------------------------------
1062   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
1063   //-----------------------------------------------------------------
1064   //
1065   if (fClustersArray==0){
1066     return;
1067   }
1068   //
1069   TParticle *particle; // pointer to a given particle
1070   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1071   //  Int_t sector,nhits;
1072   Int_t ipart;
1073   const Int_t kcmaxhits=30000;
1074   AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
1075   AliTPCFastVector & xxx = *xxxx;
1076   Int_t maxhits = kcmaxhits;
1077   //construct array for each padrow
1078   for (Int_t i=0; i<fTPCParam->GetNRow(isec);i++)
1079     fClustersArray->CreateRow(isec,i);
1080
1081   //---------------------------------------------------------------
1082   //  Get the access to the tracks
1083   //---------------------------------------------------------------
1084
1085   TTree *tH = TreeH();
1086   if (tH == 0x0)
1087    {
1088      Fatal("Hits2Clusters","Can not find TreeH in folder");
1089      return;
1090    }
1091   SetTreeAddress();
1092
1093   Stat_t ntracks = tH->GetEntries();
1094   Int_t npart = gAlice->GetMCApp()->GetNtrack();
1095   //MI change
1096   TBranch * branch=0;
1097   if (fHitType>1) branch = tH->GetBranch("TPC2");
1098   else branch = tH->GetBranch("TPC");
1099
1100   //------------------------------------------------------------
1101   // Loop over tracks
1102   //------------------------------------------------------------
1103
1104   for(Int_t track=0;track<ntracks;track++){
1105     Bool_t isInSector=kTRUE;
1106     ResetHits();
1107      isInSector = TrackInVolume(isec,track);
1108     if (!isInSector) continue;
1109     //MI change
1110     branch->GetEntry(track); // get next track
1111     //
1112     //  Get number of the TPC hits and a pointer
1113     //  to the particles
1114     // Loop over hits
1115     //
1116     Int_t currentIndex=0;
1117     Int_t lastrow=-1;  //last writen row
1118
1119     //M.I. changes
1120
1121     tpcHit = (AliTPChit*)FirstHit(-1);
1122     while(tpcHit){
1123
1124       Int_t sector=tpcHit->fSector; // sector number
1125       if(sector != isec){
1126         tpcHit = (AliTPChit*) NextHit();
1127         continue;
1128       }
1129
1130       ipart=tpcHit->Track();
1131       if (ipart<npart) particle=gAlice->GetMCApp()->Particle(ipart);
1132
1133       //find row number
1134
1135       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
1136       Int_t    index[3]={1,isec,0};
1137       Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
1138       if (currentrow<0) {tpcHit = (AliTPChit*)NextHit(); continue;}
1139       if (lastrow<0) lastrow=currentrow;
1140       if (currentrow==lastrow){
1141         if ( currentIndex>=maxhits){
1142           maxhits+=kcmaxhits;
1143           xxx.ResizeTo(4*maxhits);
1144         }
1145         xxx(currentIndex*4)=x[0];
1146         xxx(currentIndex*4+1)=x[1];
1147         xxx(currentIndex*4+2)=x[2];
1148         xxx(currentIndex*4+3)=tpcHit->fQ;
1149         currentIndex++;
1150       }
1151       else
1152         if (currentIndex>2){
1153           Float_t sumx=0;
1154           Float_t sumx2=0;
1155           Float_t sumx3=0;
1156           Float_t sumx4=0;
1157           Float_t sumy=0;
1158           Float_t sumxy=0;
1159           Float_t sumx2y=0;
1160           Float_t sumz=0;
1161           Float_t sumxz=0;
1162           Float_t sumx2z=0;
1163           Float_t sumq=0;
1164           for (Int_t index=0;index<currentIndex;index++){
1165             Float_t x,x2,x3,x4;
1166             x=x2=x3=x4=xxx(index*4);
1167             x2*=x;
1168             x3*=x2;
1169             x4*=x3;
1170             sumx+=x;
1171             sumx2+=x2;
1172             sumx3+=x3;
1173             sumx4+=x4;
1174             sumy+=xxx(index*4+1);
1175             sumxy+=xxx(index*4+1)*x;
1176             sumx2y+=xxx(index*4+1)*x2;
1177             sumz+=xxx(index*4+2);
1178             sumxz+=xxx(index*4+2)*x;
1179             sumx2z+=xxx(index*4+2)*x2;
1180             sumq+=xxx(index*4+3);
1181           }
1182           Float_t centralPad = (fTPCParam->GetNPads(isec,lastrow)-1)/2;
1183           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
1184             sumx2*(sumx*sumx3-sumx2*sumx2);
1185
1186           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
1187             sumx2*(sumxy*sumx3-sumx2y*sumx2);
1188           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
1189             sumx2*(sumxz*sumx3-sumx2z*sumx2);
1190
1191           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
1192             sumx2*(sumx*sumx2y-sumx2*sumxy);
1193           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
1194             sumx2*(sumx*sumx2z-sumx2*sumxz);
1195
1196           if (TMath::Abs(det)<0.00001){
1197              tpcHit = (AliTPChit*)NextHit();
1198             continue;
1199           }
1200
1201           Float_t y=detay/det+centralPad;
1202           Float_t z=detaz/det;
1203           Float_t by=detby/det; //y angle
1204           Float_t bz=detbz/det; //z angle
1205           sumy/=Float_t(currentIndex);
1206           sumz/=Float_t(currentIndex);
1207
1208           AliTPCClustersRow * row = (fClustersArray->GetRow(isec,lastrow));
1209           if (row!=0) {
1210             AliComplexCluster* cl = new((AliComplexCluster*)row->Append()) AliComplexCluster ;
1211             //    AliComplexCluster cl;
1212             cl->fX=z;
1213             cl->fY=y;
1214             cl->fQ=sumq;
1215             cl->fSigmaX2=bz;
1216             cl->fSigmaY2=by;
1217             cl->fTracks[0]=ipart;
1218           }
1219           currentIndex=0;
1220           lastrow=currentrow;
1221         } //end of calculating cluster for given row
1222
1223
1224       tpcHit = (AliTPChit*)NextHit();
1225     } // end of loop over hits
1226   }   // end of loop over tracks
1227   //write padrows to tree
1228   for (Int_t ii=0; ii<fTPCParam->GetNRow(isec);ii++) {
1229     fClustersArray->StoreRow(isec,ii);
1230     fClustersArray->ClearRow(isec,ii);
1231   }
1232   xxxx->Delete();
1233
1234 }
1235
1236
1237
1238 //______________________________________________________________________
1239 AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
1240 {
1241   return new AliTPCDigitizer(manager);
1242 }
1243 //__
1244 void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)
1245 {
1246   //create digits from summable digits
1247   GenerNoise(500000); //create teble with noise
1248
1249   //conect tree with sSDigits
1250   TTree *t = fLoader->TreeS();
1251
1252   if (t == 0x0)
1253    {
1254      fLoader->LoadSDigits("READ");
1255      t = fLoader->TreeS();
1256      if (t == 0x0)
1257       {
1258         Error("SDigits2Digits2","Can not get input TreeS");
1259         return;
1260       }
1261    }
1262
1263   if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D");
1264
1265   AliSimDigits digarr, *dummy=&digarr;
1266   TBranch* sdb = t->GetBranch("Segment");
1267   if (sdb == 0x0)
1268    {
1269      Error("SDigits2Digits2","Can not find branch with segments in TreeS.");
1270      return;
1271    }
1272
1273   sdb->SetAddress(&dummy);
1274
1275   Stat_t nentries = t->GetEntries();
1276
1277   // set zero suppression
1278
1279   fTPCParam->SetZeroSup(2);
1280
1281   // get zero suppression
1282
1283   Int_t zerosup = fTPCParam->GetZeroSup();
1284
1285   //make tree with digits
1286
1287   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1288   arr->SetClass("AliSimDigits");
1289   arr->Setup(fTPCParam);
1290   arr->MakeTree(fLoader->TreeD());
1291
1292   AliTPCParam * par = fTPCParam;
1293
1294   //Loop over segments of the TPC
1295
1296   for (Int_t n=0; n<nentries; n++) {
1297     t->GetEvent(n);
1298     Int_t sec, row;
1299     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
1300       cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
1301       continue;
1302     }
1303     if (!IsSectorActive(sec))
1304      {
1305 //       cout<<n<<" NOT Active \n";
1306        continue;
1307      }
1308     else
1309      {
1310 //       cout<<n<<" Active \n";
1311      }
1312     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
1313     Int_t nrows = digrow->GetNRows();
1314     Int_t ncols = digrow->GetNCols();
1315
1316     digrow->ExpandBuffer();
1317     digarr.ExpandBuffer();
1318     digrow->ExpandTrackBuffer();
1319     digarr.ExpandTrackBuffer();
1320
1321
1322     Short_t * pamp0 = digarr.GetDigits();
1323     Int_t   * ptracks0 = digarr.GetTracks();
1324     Short_t * pamp1 = digrow->GetDigits();
1325     Int_t   * ptracks1 = digrow->GetTracks();
1326     Int_t  nelems =nrows*ncols;
1327     Int_t saturation = fTPCParam->GetADCSat();
1328     //use internal structure of the AliDigits - for speed reason
1329     //if you cahnge implementation
1330     //of the Alidigits - it must be rewriten -
1331     for (Int_t i= 0; i<nelems; i++){
1332       //      Float_t q = *pamp0;
1333       //q/=16.;  //conversion faktor
1334       //Float_t noise= GetNoise();
1335       //q+=noise;
1336       //q= TMath::Nint(q);
1337       Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
1338       if (q>zerosup){
1339         if (q>saturation) q=saturation;
1340         *pamp1=(Short_t)q;
1341         //if (ptracks0[0]==0)
1342         //  ptracks1[0]=1;
1343         //else
1344         ptracks1[0]=ptracks0[0];
1345         ptracks1[nelems]=ptracks0[nelems];
1346         ptracks1[2*nelems]=ptracks0[2*nelems];
1347       }
1348       pamp0++;
1349       pamp1++;
1350       ptracks0++;
1351       ptracks1++;
1352     }
1353
1354     arr->StoreRow(sec,row);
1355     arr->ClearRow(sec,row);
1356     // cerr<<sec<<"\t"<<row<<"\n";
1357   }
1358
1359
1360   //write results
1361   fLoader->WriteDigits("OVERWRITE");
1362
1363   delete arr;
1364 }
1365 //__________________________________________________________________
1366 void AliTPC::SetDefaults(){
1367   //
1368   // setting the defaults
1369   //
1370
1371    cerr<<"Setting default parameters...\n";
1372
1373   // Set response functions
1374
1375   //
1376   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
1377   rl->CdGAFile();
1378   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1379   if(param){
1380     printf("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...\n");
1381     delete param;
1382     param = new AliTPCParamSR();
1383   }
1384   else {
1385     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1386   }
1387   if(!param){
1388     printf("No TPC parameters found\n");
1389     exit(4);
1390   }
1391
1392
1393   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1394   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1395   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;
1396   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1397   rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1398   rf->SetOffset(3*param->GetZSigma());
1399   rf->Update();
1400
1401   TDirectory *savedir=gDirectory;
1402   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1403   if (!f->IsOpen()) {
1404     cerr<<"Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !\n" ;
1405      exit(3);
1406   }
1407
1408   TString s;
1409   prfinner->Read("prf_07504_Gati_056068_d02");
1410   //PH Set different names
1411   s=prfinner->GetGRF()->GetName();
1412   s+="in";
1413   prfinner->GetGRF()->SetName(s.Data());
1414
1415   prfouter1->Read("prf_10006_Gati_047051_d03");
1416   s=prfouter1->GetGRF()->GetName();
1417   s+="out1";
1418   prfouter1->GetGRF()->SetName(s.Data());
1419
1420   prfouter2->Read("prf_15006_Gati_047051_d03");
1421   s=prfouter2->GetGRF()->GetName();
1422   s+="out2";
1423   prfouter2->GetGRF()->SetName(s.Data());
1424
1425   f->Close();
1426   savedir->cd();
1427
1428   param->SetInnerPRF(prfinner);
1429   param->SetOuter1PRF(prfouter1);
1430   param->SetOuter2PRF(prfouter2);
1431   param->SetTimeRF(rf);
1432
1433   // set fTPCParam
1434
1435   SetParam(param);
1436
1437
1438   fDefaults = 1;
1439
1440 }
1441 //__________________________________________________________________
1442 void AliTPC::Hits2Digits()
1443 {
1444   //
1445   // creates digits from hits
1446   //
1447
1448   fLoader->LoadHits("read");
1449   fLoader->LoadDigits("recreate");
1450   AliRunLoader* runLoader = fLoader->GetRunLoader();
1451
1452   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1453     runLoader->GetEvent(iEvent);
1454     SetActiveSectors();
1455     Hits2Digits(iEvent);
1456   }
1457
1458   fLoader->UnloadHits();
1459   fLoader->UnloadDigits();
1460 }
1461 //__________________________________________________________________
1462 void AliTPC::Hits2Digits(Int_t eventnumber)
1463 {
1464  //----------------------------------------------------
1465  // Loop over all sectors for a single event
1466  //----------------------------------------------------
1467   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
1468   rl->GetEvent(eventnumber);
1469   if (fLoader->TreeH() == 0x0)
1470    {
1471      if(fLoader->LoadHits())
1472       {
1473         Error("Hits2Digits","Can not load hits.");
1474       }
1475    }
1476   SetTreeAddress();
1477
1478   if (fLoader->TreeD() == 0x0 )
1479    {
1480      fLoader->MakeTree("D");
1481      if (fLoader->TreeD() == 0x0 )
1482       {
1483        Error("Hits2Digits","Can not get TreeD");
1484        return;
1485       }
1486    }
1487
1488   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1489   GenerNoise(500000); //create teble with noise
1490
1491   //setup TPCDigitsArray
1492
1493   if(GetDigitsArray()) delete GetDigitsArray();
1494
1495   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1496   arr->SetClass("AliSimDigits");
1497   arr->Setup(fTPCParam);
1498
1499   arr->MakeTree(fLoader->TreeD());
1500   SetDigitsArray(arr);
1501
1502   fDigitsSwitch=0; // standard digits
1503
1504   cerr<<"Digitizing TPC -- normal digits...\n";
1505
1506  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1507   if (IsSectorActive(isec))
1508    {
1509     if (fDebug) Info("Hits2Digits","Sector %d is active.",isec);
1510     Hits2DigitsSector(isec);
1511    }
1512   else
1513    {
1514     if (fDebug) Info("Hits2Digits","Sector %d is NOT active.",isec);
1515    }
1516
1517   fLoader->WriteDigits("OVERWRITE");
1518
1519 //this line prevents the crash in the similar one
1520 //on the beginning of this method
1521 //destructor attempts to reset the tree, which is deleted by the loader
1522 //need to be redesign
1523  if(GetDigitsArray()) delete GetDigitsArray();
1524  SetDigitsArray(0x0);
1525
1526 }
1527
1528 //__________________________________________________________________
1529 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1530 {
1531
1532   //-----------------------------------------------------------
1533   //   summable digits - 16 bit "ADC", no noise, no saturation
1534   //-----------------------------------------------------------
1535
1536  //----------------------------------------------------
1537  // Loop over all sectors for a single event
1538  //----------------------------------------------------
1539 //  AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::fgkRunLoaderName);
1540
1541   AliRunLoader* rl = fLoader->GetRunLoader();
1542
1543   rl->GetEvent(eventnumber);
1544   if (fLoader->TreeH() == 0x0)
1545    {
1546      if(fLoader->LoadHits())
1547       {
1548         Error("Hits2Digits","Can not load hits.");
1549         return;
1550       }
1551    }
1552   SetTreeAddress();
1553
1554
1555   if (fLoader->TreeS() == 0x0 )
1556    {
1557      fLoader->MakeTree("S");
1558    }
1559
1560   if(fDefaults == 0) SetDefaults();
1561
1562   GenerNoise(500000); //create table with noise
1563   //setup TPCDigitsArray
1564
1565   if(GetDigitsArray()) delete GetDigitsArray();
1566
1567
1568   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1569   arr->SetClass("AliSimDigits");
1570   arr->Setup(fTPCParam);
1571   arr->MakeTree(fLoader->TreeS());
1572
1573   SetDigitsArray(arr);
1574
1575   cerr<<"Digitizing TPC -- summable digits...\n";
1576
1577   fDigitsSwitch=1; // summable digits
1578
1579     // set zero suppression to "0"
1580
1581   fTPCParam->SetZeroSup(0);
1582
1583  for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1584   if (IsSectorActive(isec))
1585    {
1586 //    cout<<"Sector "<<isec<<" is active\n";
1587     Hits2DigitsSector(isec);
1588    }
1589
1590  fLoader->WriteSDigits("OVERWRITE");
1591
1592 //this line prevents the crash in the similar one
1593 //on the beginning of this method
1594 //destructor attempts to reset the tree, which is deleted by the loader
1595 //need to be redesign
1596  if(GetDigitsArray()) delete GetDigitsArray();
1597  SetDigitsArray(0x0);
1598 }
1599 //__________________________________________________________________
1600
1601 void AliTPC::Hits2SDigits()
1602 {
1603
1604   //-----------------------------------------------------------
1605   //   summable digits - 16 bit "ADC", no noise, no saturation
1606   //-----------------------------------------------------------
1607
1608   fLoader->LoadHits("read");
1609   fLoader->LoadSDigits("recreate");
1610   AliRunLoader* runLoader = fLoader->GetRunLoader();
1611
1612   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1613     runLoader->GetEvent(iEvent);
1614     SetTreeAddress();
1615     SetActiveSectors();
1616     Hits2SDigits2(iEvent);
1617   }
1618
1619   fLoader->UnloadHits();
1620   fLoader->UnloadSDigits();
1621 }
1622 //_____________________________________________________________________________
1623
1624 void AliTPC::Hits2DigitsSector(Int_t isec)
1625 {
1626   //-------------------------------------------------------------------
1627   // TPC conversion from hits to digits.
1628   //-------------------------------------------------------------------
1629
1630   //-----------------------------------------------------------------
1631   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1632   //-----------------------------------------------------------------
1633
1634   //-------------------------------------------------------
1635   //  Get the access to the track hits
1636   //-------------------------------------------------------
1637
1638   // check if the parameters are set - important if one calls this method
1639   // directly, not from the Hits2Digits
1640
1641   if(fDefaults == 0) SetDefaults();
1642
1643   TTree *tH = TreeH(); // pointer to the hits tree
1644   if (tH == 0x0)
1645    {
1646      Fatal("Hits2DigitsSector","Can not find TreeH in folder");
1647      return;
1648    }
1649
1650   Stat_t ntracks = tH->GetEntries();
1651
1652   if( ntracks > 0){
1653
1654   //-------------------------------------------
1655   //  Only if there are any tracks...
1656   //-------------------------------------------
1657
1658     TObjArray **row;
1659
1660     //printf("*** Processing sector number %d ***\n",isec);
1661
1662       Int_t nrows =fTPCParam->GetNRow(isec);
1663
1664       row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1665
1666       MakeSector(isec,nrows,tH,ntracks,row);
1667
1668       //--------------------------------------------------------
1669       //   Digitize this sector, row by row
1670       //   row[i] is the pointer to the TObjArray of AliTPCFastVectors,
1671       //   each one containing electrons accepted on this
1672       //   row, assigned into tracks
1673       //--------------------------------------------------------
1674
1675       Int_t i;
1676
1677       if (fDigitsArray->GetTree()==0)
1678        {
1679          Fatal("Hits2DigitsSector","Tree not set in fDigitsArray");
1680        }
1681
1682       for (i=0;i<nrows;i++){
1683
1684         AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1685
1686         DigitizeRow(i,isec,row);
1687
1688         fDigitsArray->StoreRow(isec,i);
1689
1690         Int_t ndig = dig->GetDigitSize();
1691
1692         if (gDebug > 10)
1693         printf("*** Sector, row, compressed digits %d %d %d ***\n",isec,i,ndig);
1694
1695         fDigitsArray->ClearRow(isec,i);
1696
1697
1698        } // end of the sector digitization
1699
1700       for(i=0;i<nrows+2;i++){
1701         row[i]->Delete();
1702         delete row[i];
1703       }
1704
1705        delete [] row; // delete the array of pointers to TObjArray-s
1706
1707   } // ntracks >0
1708
1709 } // end of Hits2DigitsSector
1710
1711
1712 //_____________________________________________________________________________
1713 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1714 {
1715   //-----------------------------------------------------------
1716   // Single row digitization, coupling from the neighbouring
1717   // rows taken into account
1718   //-----------------------------------------------------------
1719
1720   //-----------------------------------------------------------------
1721   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1722   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1723   //-----------------------------------------------------------------
1724
1725
1726   Float_t zerosup = fTPCParam->GetZeroSup();
1727   //  Int_t nrows =fTPCParam->GetNRow(isec);
1728   fCurrentIndex[1]= isec;
1729
1730
1731   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1732   Int_t nofTbins = fTPCParam->GetMaxTBin();
1733   Int_t indexRange[4];
1734   //
1735   //  Integrated signal for this row
1736   //  and a single track signal
1737   //
1738
1739   AliTPCFastMatrix *m1 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // integrated
1740   AliTPCFastMatrix *m2 = new AliTPCFastMatrix(0,nofPads,0,nofTbins); // single
1741   //
1742   AliTPCFastMatrix &total  = *m1;
1743
1744   //  Array of pointers to the label-signal list
1745
1746   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1747   Float_t  **pList = new Float_t* [nofDigits];
1748
1749   Int_t lp;
1750   Int_t i1;
1751   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1752   //
1753   //calculate signal
1754   //
1755   //Int_t row1 = TMath::Max(irow-fTPCParam->GetNCrossRows(),0);
1756   //Int_t row2 = TMath::Min(irow+fTPCParam->GetNCrossRows(),nrows-1);
1757   Int_t row1=irow;
1758   Int_t row2=irow+2;
1759   for (Int_t row= row1;row<=row2;row++){
1760     Int_t nTracks= rows[row]->GetEntries();
1761     for (i1=0;i1<nTracks;i1++){
1762       fCurrentIndex[2]= row;
1763       fCurrentIndex[3]=irow+1;
1764       if (row==irow+1){
1765         m2->Zero();  // clear single track signal matrix
1766         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1767         GetList(trackLabel,nofPads,m2,indexRange,pList);
1768       }
1769       else   GetSignal(rows[row],i1,0,m1,indexRange);
1770     }
1771   }
1772
1773   Int_t tracks[3];
1774
1775   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1776   Int_t gi=-1;
1777   Float_t fzerosup = zerosup+0.5;
1778   for(Int_t it=0;it<nofTbins;it++){
1779     Float_t *pq = &(total.UncheckedAt(0,it));
1780     for(Int_t ip=0;ip<nofPads;ip++){
1781       gi++;
1782       Float_t q=*pq;
1783       pq++;
1784       if(fDigitsSwitch == 0){
1785         q+=GetNoise();
1786         if(q <=fzerosup) continue; // do not fill zeros
1787         q = TMath::Nint(q);
1788         if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();  // saturation
1789
1790       }
1791
1792       else {
1793        if(q <= 0.) continue; // do not fill zeros
1794        if(q>2000.) q=2000.;
1795        q *= 16.;
1796        q = TMath::Nint(q);
1797       }
1798
1799       //
1800       //  "real" signal or electronic noise (list = -1)?
1801       //
1802
1803       for(Int_t j1=0;j1<3;j1++){
1804         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1805       }
1806
1807 //Begin_Html
1808 /*
1809   <A NAME="AliDigits"></A>
1810   using of AliDigits object
1811 */
1812 //End_Html
1813       dig->SetDigitFast((Short_t)q,it,ip);
1814       if (fDigitsArray->IsSimulated())
1815         {
1816          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1817          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1818          ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1819         }
1820
1821
1822     } // end of loop over time buckets
1823   }  // end of lop over pads
1824
1825   //
1826   //  This row has been digitized, delete nonused stuff
1827   //
1828
1829   for(lp=0;lp<nofDigits;lp++){
1830     if(pList[lp]) delete [] pList[lp];
1831   }
1832
1833   delete [] pList;
1834
1835   delete m1;
1836   delete m2;
1837   //  delete m3;
1838
1839 } // end of DigitizeRow
1840
1841 //_____________________________________________________________________________
1842
1843 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
1844              AliTPCFastMatrix *m1, AliTPCFastMatrix *m2,Int_t *indexRange)
1845 {
1846
1847   //---------------------------------------------------------------
1848   //  Calculates 2-D signal (pad,time) for a single track,
1849   //  returns a pointer to the signal matrix and the track label
1850   //  No digitization is performed at this level!!!
1851   //---------------------------------------------------------------
1852
1853   //-----------------------------------------------------------------
1854   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1855   // Modified: Marian Ivanov
1856   //-----------------------------------------------------------------
1857
1858   AliTPCFastVector *tv;
1859
1860   tv = (AliTPCFastVector*)p1->At(ntr); // pointer to a track
1861   AliTPCFastVector &v = *tv;
1862
1863   Float_t label = v(0);
1864   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1)-1)/2;
1865
1866   Int_t nElectrons = (tv->GetNrows()-1)/4;
1867   indexRange[0]=9999; // min pad
1868   indexRange[1]=-1; // max pad
1869   indexRange[2]=9999; //min time
1870   indexRange[3]=-1; // max time
1871
1872   AliTPCFastMatrix &signal = *m1;
1873   AliTPCFastMatrix &total = *m2;
1874   //
1875   //  Loop over all electrons
1876   //
1877   for(Int_t nel=0; nel<nElectrons; nel++){
1878     Int_t idx=nel*4;
1879     Float_t aval =  v(idx+4);
1880     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1881     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1882     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
1883
1884     Int_t *index = fTPCParam->GetResBin(0);
1885     Float_t *weight = & (fTPCParam->GetResWeight(0));
1886
1887     if (n>0) for (Int_t i =0; i<n; i++){
1888        Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1889
1890          if (pad>=0){
1891          Int_t time=index[2];
1892          Float_t qweight = *(weight)*eltoadcfac;
1893
1894          if (m1!=0) signal.UncheckedAt(pad,time)+=qweight;
1895          total.UncheckedAt(pad,time)+=qweight;
1896          if (indexRange[0]>pad) indexRange[0]=pad;
1897          if (indexRange[1]<pad) indexRange[1]=pad;
1898          if (indexRange[2]>time) indexRange[2]=time;
1899          if (indexRange[3]<time) indexRange[3]=time;
1900
1901          index+=3;
1902          weight++;
1903
1904        }
1905     }
1906   } // end of loop over electrons
1907
1908   return label; // returns track label when finished
1909 }
1910
1911 //_____________________________________________________________________________
1912 void AliTPC::GetList(Float_t label,Int_t np,AliTPCFastMatrix *m,
1913                      Int_t *indexRange, Float_t **pList)
1914 {
1915   //----------------------------------------------------------------------
1916   //  Updates the list of tracks contributing to digits for a given row
1917   //----------------------------------------------------------------------
1918
1919   //-----------------------------------------------------------------
1920   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1921   //-----------------------------------------------------------------
1922
1923   AliTPCFastMatrix &signal = *m;
1924
1925   // lop over nonzero digits
1926
1927   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1928     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1929
1930
1931         // accept only the contribution larger than 500 electrons (1/2 s_noise)
1932
1933         if(signal(ip,it)<0.5) continue;
1934
1935
1936         Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1937
1938         if(!pList[globalIndex]){
1939
1940           //
1941           // Create new list (6 elements - 3 signals and 3 labels),
1942           //
1943
1944           pList[globalIndex] = new Float_t [6];
1945
1946           // set list to -1
1947
1948           *pList[globalIndex] = -1.;
1949           *(pList[globalIndex]+1) = -1.;
1950           *(pList[globalIndex]+2) = -1.;
1951           *(pList[globalIndex]+3) = -1.;
1952           *(pList[globalIndex]+4) = -1.;
1953           *(pList[globalIndex]+5) = -1.;
1954
1955
1956           *pList[globalIndex] = label;
1957           *(pList[globalIndex]+3) = signal(ip,it);
1958         }
1959         else{
1960
1961           // check the signal magnitude
1962
1963           Float_t highest = *(pList[globalIndex]+3);
1964           Float_t middle = *(pList[globalIndex]+4);
1965           Float_t lowest = *(pList[globalIndex]+5);
1966
1967           //
1968           //  compare the new signal with already existing list
1969           //
1970
1971           if(signal(ip,it)<lowest) continue; // neglect this track
1972
1973           //
1974
1975           if (signal(ip,it)>highest){
1976             *(pList[globalIndex]+5) = middle;
1977             *(pList[globalIndex]+4) = highest;
1978             *(pList[globalIndex]+3) = signal(ip,it);
1979
1980             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1981             *(pList[globalIndex]+1) = *pList[globalIndex];
1982             *pList[globalIndex] = label;
1983           }
1984           else if (signal(ip,it)>middle){
1985             *(pList[globalIndex]+5) = middle;
1986             *(pList[globalIndex]+4) = signal(ip,it);
1987
1988             *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1989             *(pList[globalIndex]+1) = label;
1990           }
1991           else{
1992             *(pList[globalIndex]+5) = signal(ip,it);
1993             *(pList[globalIndex]+2) = label;
1994           }
1995         }
1996
1997     } // end of loop over pads
1998   } // end of loop over time bins
1999
2000
2001
2002 }//end of GetList
2003 //___________________________________________________________________
2004 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
2005                         Stat_t ntracks,TObjArray **row)
2006 {
2007
2008   //-----------------------------------------------------------------
2009   // Prepares the sector digitization, creates the vectors of
2010   // tracks for each row of this sector. The track vector
2011   // contains the track label and the position of electrons.
2012   //-----------------------------------------------------------------
2013
2014   //-----------------------------------------------------------------
2015   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2016   //-----------------------------------------------------------------
2017
2018   Float_t gasgain = fTPCParam->GetGasGain();
2019   Int_t i;
2020   Float_t xyz[4];
2021
2022   AliTPChit *tpcHit; // pointer to a sigle TPC hit
2023   //MI change
2024   TBranch * branch=0;
2025   if (fHitType>1) branch = TH->GetBranch("TPC2");
2026   else branch = TH->GetBranch("TPC");
2027
2028
2029   //----------------------------------------------
2030   // Create TObjArray-s, one for each row,
2031   // each TObjArray will store the AliTPCFastVectors
2032   // of electrons, one AliTPCFastVectors per each track.
2033   //----------------------------------------------
2034
2035   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
2036   AliTPCFastVector **tracks = new AliTPCFastVector* [nrows+2]; //pointers to the track vectors
2037
2038   for(i=0; i<nrows+2; i++){
2039     row[i] = new TObjArray;
2040     nofElectrons[i]=0;
2041     tracks[i]=0;
2042   }
2043
2044
2045
2046   //--------------------------------------------------------------------
2047   //  Loop over tracks, the "track" contains the full history
2048   //--------------------------------------------------------------------
2049
2050   Int_t previousTrack,currentTrack;
2051   previousTrack = -1; // nothing to store so far!
2052
2053   for(Int_t track=0;track<ntracks;track++){
2054     Bool_t isInSector=kTRUE;
2055     ResetHits();
2056     isInSector = TrackInVolume(isec,track);
2057     if (!isInSector) continue;
2058     //MI change
2059     branch->GetEntry(track); // get next track
2060
2061     //M.I. changes
2062
2063     tpcHit = (AliTPChit*)FirstHit(-1);
2064
2065     //--------------------------------------------------------------
2066     //  Loop over hits
2067     //--------------------------------------------------------------
2068
2069
2070     while(tpcHit){
2071
2072       Int_t sector=tpcHit->fSector; // sector number
2073       if(sector != isec){
2074         tpcHit = (AliTPChit*) NextHit();
2075         continue;
2076       }
2077
2078         currentTrack = tpcHit->Track(); // track number
2079
2080
2081         if(currentTrack != previousTrack){
2082
2083            // store already filled fTrack
2084
2085            for(i=0;i<nrows+2;i++){
2086              if(previousTrack != -1){
2087                if(nofElectrons[i]>0){
2088                  AliTPCFastVector &v = *tracks[i];
2089                  v(0) = previousTrack;
2090                  tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2091                  row[i]->Add(tracks[i]);
2092                }
2093                else{
2094                  delete tracks[i]; // delete empty AliTPCFastVector
2095                  tracks[i]=0;
2096                }
2097              }
2098
2099              nofElectrons[i]=0;
2100              tracks[i] = new AliTPCFastVector(481); // AliTPCFastVectors for the next fTrack
2101
2102            } // end of loop over rows
2103
2104            previousTrack=currentTrack; // update track label
2105         }
2106
2107         Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2108
2109        //---------------------------------------------------
2110        //  Calculate the electron attachment probability
2111        //---------------------------------------------------
2112
2113
2114         Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
2115                                                         /fTPCParam->GetDriftV();
2116         // in microseconds!
2117         Float_t attProb = fTPCParam->GetAttCoef()*
2118           fTPCParam->GetOxyCont()*time; //  fraction!
2119
2120         //-----------------------------------------------
2121         //  Loop over electrons
2122         //-----------------------------------------------
2123         Int_t index[3];
2124         index[1]=isec;
2125         for(Int_t nel=0;nel<qI;nel++){
2126           // skip if electron lost due to the attachment
2127           if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
2128           xyz[0]=tpcHit->X();
2129           xyz[1]=tpcHit->Y();
2130           xyz[2]=tpcHit->Z();
2131           //
2132           // protection for the nonphysical avalanche size (10**6 maximum)
2133           //
2134           Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2135           xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
2136           index[0]=1;
2137
2138           TransportElectron(xyz,index);
2139           Int_t rowNumber;
2140           fTPCParam->GetPadRow(xyz,index);
2141           // row 0 - cross talk from the innermost row
2142           // row fNRow+1 cross talk from the outermost row
2143           rowNumber = index[2]+1;
2144           //transform position to local digit coordinates
2145           //relative to nearest pad row
2146           if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
2147           Float_t x1,y1;
2148           if (isec <fTPCParam->GetNInnerSector()) {
2149             x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
2150             y1 = fTPCParam->GetYInner(rowNumber);
2151           }
2152           else{
2153             x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
2154             y1 = fTPCParam->GetYOuter(rowNumber);
2155           }
2156           // gain inefficiency at the wires edges - linear
2157           x1=TMath::Abs(x1);
2158           y1-=1.;
2159           if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));
2160
2161           nofElectrons[rowNumber]++;
2162           //----------------------------------
2163           // Expand vector if necessary
2164           //----------------------------------
2165           if(nofElectrons[rowNumber]>120){
2166             Int_t range = tracks[rowNumber]->GetNrows();
2167             if((nofElectrons[rowNumber])>(range-1)/4){
2168
2169               tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
2170             }
2171           }
2172
2173           AliTPCFastVector &v = *tracks[rowNumber];
2174           Int_t idx = 4*nofElectrons[rowNumber]-3;
2175           Real_t * position = &(((AliTPCFastVector&)v).UncheckedAt(idx)); //make code faster
2176           memcpy(position,xyz,4*sizeof(Float_t));
2177
2178         } // end of loop over electrons
2179
2180         tpcHit = (AliTPChit*)NextHit();
2181
2182       } // end of loop over hits
2183     } // end of loop over tracks
2184
2185     //
2186     //   store remaining track (the last one) if not empty
2187     //
2188
2189      for(i=0;i<nrows+2;i++){
2190        if(nofElectrons[i]>0){
2191           AliTPCFastVector &v = *tracks[i];
2192           v(0) = previousTrack;
2193           tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
2194           row[i]->Add(tracks[i]);
2195         }
2196         else{
2197           delete tracks[i];
2198           tracks[i]=0;
2199         }
2200       }
2201
2202           delete [] tracks;
2203           delete [] nofElectrons;
2204
2205
2206 } // end of MakeSector
2207
2208
2209 //_____________________________________________________________________________
2210 void AliTPC::Init()
2211 {
2212   //
2213   // Initialise TPC detector after definition of geometry
2214   //
2215   Int_t i;
2216   //
2217   if(fDebug) {
2218     printf("\n%s: ",ClassName());
2219     for(i=0;i<35;i++) printf("*");
2220     printf(" TPC_INIT ");
2221     for(i=0;i<35;i++) printf("*");
2222     printf("\n%s: ",ClassName());
2223     //
2224     for(i=0;i<80;i++) printf("*");
2225     printf("\n");
2226   }
2227 }
2228
2229 //_____________________________________________________________________________
2230 void AliTPC::MakeBranch(Option_t* option)
2231 {
2232   //
2233   // Create Tree branches for the TPC.
2234   //
2235   if(GetDebug()) Info("MakeBranch","");
2236   Int_t buffersize = 4000;
2237   char branchname[10];
2238   sprintf(branchname,"%s",GetName());
2239
2240   const char *h = strstr(option,"H");
2241
2242   if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2243
2244   AliDetector::MakeBranch(option);
2245
2246   const char *d = strstr(option,"D");
2247
2248   if (fDigits   && fLoader->TreeD() && d)
2249    {
2250       MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0);
2251    }
2252
2253   if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000
2254 }
2255
2256 //_____________________________________________________________________________
2257 void AliTPC::ResetDigits()
2258 {
2259   //
2260   // Reset number of digits and the digits array for this detector
2261   //
2262   fNdigits   = 0;
2263   if (fDigits)   fDigits->Clear();
2264 }
2265
2266 //_____________________________________________________________________________
2267 void AliTPC::SetSecAL(Int_t sec)
2268 {
2269   //---------------------------------------------------
2270   // Activate/deactivate selection for lower sectors
2271   //---------------------------------------------------
2272
2273   //-----------------------------------------------------------------
2274   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2275   //-----------------------------------------------------------------
2276   fSecAL = sec;
2277 }
2278
2279 //_____________________________________________________________________________
2280 void AliTPC::SetSecAU(Int_t sec)
2281 {
2282   //----------------------------------------------------
2283   // Activate/deactivate selection for upper sectors
2284   //---------------------------------------------------
2285
2286   //-----------------------------------------------------------------
2287   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2288   //-----------------------------------------------------------------
2289   fSecAU = sec;
2290 }
2291
2292 //_____________________________________________________________________________
2293 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
2294 {
2295   //----------------------------------------
2296   // Select active lower sectors
2297   //----------------------------------------
2298
2299   //-----------------------------------------------------------------
2300   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2301   //-----------------------------------------------------------------
2302
2303   fSecLows[0] = s1;
2304   fSecLows[1] = s2;
2305   fSecLows[2] = s3;
2306   fSecLows[3] = s4;
2307   fSecLows[4] = s5;
2308   fSecLows[5] = s6;
2309 }
2310
2311 //_____________________________________________________________________________
2312 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
2313                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
2314                        Int_t s11 , Int_t s12)
2315 {
2316   //--------------------------------
2317   // Select active upper sectors
2318   //--------------------------------
2319
2320   //-----------------------------------------------------------------
2321   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2322   //-----------------------------------------------------------------
2323
2324   fSecUps[0] = s1;
2325   fSecUps[1] = s2;
2326   fSecUps[2] = s3;
2327   fSecUps[3] = s4;
2328   fSecUps[4] = s5;
2329   fSecUps[5] = s6;
2330   fSecUps[6] = s7;
2331   fSecUps[7] = s8;
2332   fSecUps[8] = s9;
2333   fSecUps[9] = s10;
2334   fSecUps[10] = s11;
2335   fSecUps[11] = s12;
2336 }
2337
2338 //_____________________________________________________________________________
2339 void AliTPC::SetSens(Int_t sens)
2340 {
2341
2342   //-------------------------------------------------------------
2343   // Activates/deactivates the sensitive strips at the center of
2344   // the pad row -- this is for the space-point resolution calculations
2345   //-------------------------------------------------------------
2346
2347   //-----------------------------------------------------------------
2348   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2349   //-----------------------------------------------------------------
2350
2351   fSens = sens;
2352 }
2353
2354
2355 void AliTPC::SetSide(Float_t side=0.)
2356 {
2357   // choice of the TPC side
2358
2359   fSide = side;
2360
2361 }
2362 //____________________________________________________________________________
2363 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
2364                            Float_t p2,Float_t p3)
2365 {
2366
2367   // gax mixture definition
2368
2369  fNoComp = nc;
2370
2371  fMixtComp[0]=c1;
2372  fMixtComp[1]=c2;
2373  fMixtComp[2]=c3;
2374
2375  fMixtProp[0]=p1;
2376  fMixtProp[1]=p2;
2377  fMixtProp[2]=p3;
2378
2379
2380 }
2381 //_____________________________________________________________________________
2382
2383 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2384 {
2385   //
2386   // electron transport taking into account:
2387   // 1. diffusion,
2388   // 2.ExB at the wires
2389   // 3. nonisochronity
2390   //
2391   // xyz and index must be already transformed to system 1
2392   //
2393
2394   fTPCParam->Transform1to2(xyz,index);
2395
2396   //add diffusion
2397   Float_t driftl=xyz[2];
2398   if(driftl<0.01) driftl=0.01;
2399   driftl=TMath::Sqrt(driftl);
2400   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2401   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2402   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2403   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2404   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2405
2406   // ExB
2407
2408   if (fTPCParam->GetMWPCReadout()==kTRUE){
2409     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2410     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2411   }
2412   //add nonisochronity (not implemented yet)
2413 }
2414
2415 ClassImp(AliTPCdigit)
2416
2417 //_____________________________________________________________________________
2418 AliTPCdigit::AliTPCdigit(Int_t *tracks, Int_t *digits):
2419   AliDigit(tracks)
2420 {
2421   //
2422   // Creates a TPC digit object
2423   //
2424   fSector     = digits[0];
2425   fPadRow     = digits[1];
2426   fPad        = digits[2];
2427   fTime       = digits[3];
2428   fSignal     = digits[4];
2429 }
2430
2431
2432 ClassImp(AliTPChit)
2433
2434 //_____________________________________________________________________________
2435 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2436 AliHit(shunt,track)
2437 {
2438   //
2439   // Creates a TPC hit object
2440   //
2441   fSector     = vol[0];
2442   fPadRow     = vol[1];
2443   fX          = hits[0];
2444   fY          = hits[1];
2445   fZ          = hits[2];
2446   fQ          = hits[3];
2447 }
2448
2449 //________________________________________________________________________
2450 // Additional code because of the AliTPCTrackHitsV2
2451
2452 void AliTPC::MakeBranch2(Option_t *option,const char */*file*/)
2453 {
2454   //
2455   // Create a new branch in the current Root Tree
2456   // The branch of fHits is automatically split
2457   // MI change 14.09.2000
2458   if(GetDebug()) Info("MakeBranch2","");
2459   if (fHitType<2) return;
2460   char branchname[10];
2461   sprintf(branchname,"%s2",GetName());
2462   //
2463   // Get the pointer to the header
2464   const char *cH = strstr(option,"H");
2465   //
2466   if (fTrackHits   && TreeH() && cH && fHitType&4)
2467    {
2468     if(GetDebug()) Info("MakeBranch2","Making branch for Type 4 Hits");
2469     TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2470    }
2471
2472   if (fTrackHitsOld   && TreeH() && cH && fHitType&2)
2473    {
2474     if(GetDebug()) Info("MakeBranch2","Making branch for Type 2 Hits");
2475     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2476                                                    TreeH(),fBufferSize,99);
2477     TreeH()->GetListOfBranches()->Add(branch);
2478    }
2479 }
2480
2481 void AliTPC::SetTreeAddress()
2482 {
2483 //Sets tree address for hits
2484   if (fHitType<=1)
2485    {
2486      if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2487      AliDetector::SetTreeAddress();
2488    }
2489   if (fHitType>1) SetTreeAddress2();
2490 }
2491
2492 void AliTPC::SetTreeAddress2()
2493 {
2494   //
2495   // Set branch address for the TrackHits Tree
2496   //
2497   if(GetDebug()) Info("SetTreeAddress2","");
2498
2499   TBranch *branch;
2500   char branchname[20];
2501   sprintf(branchname,"%s2",GetName());
2502   //
2503   // Branch address for hit tree
2504   TTree *treeH = TreeH();
2505   if ((treeH)&&(fHitType&4)) {
2506     branch = treeH->GetBranch(branchname);
2507     if (branch)
2508      {
2509        branch->SetAddress(&fTrackHits);
2510        if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Setting");
2511      }
2512     else
2513     if (GetDebug()) Info("SetTreeAddress2","fHitType&4 Failed (can not find branch)");
2514
2515   }
2516   if ((treeH)&&(fHitType&2)) {
2517     branch = treeH->GetBranch(branchname);
2518     if (branch)
2519      {
2520        branch->SetAddress(&fTrackHitsOld);
2521        if (GetDebug()) Info("SetTreeAddress2","fHitType&2 Setting");
2522      }
2523     else if (GetDebug())
2524       Info("SetTreeAddress2","fHitType&2 Failed (can not find branch)");
2525   }
2526   //set address to TREETR
2527
2528   TTree *treeTR = TreeTR();
2529   if (treeTR && fTrackReferences) {
2530     branch = treeTR->GetBranch(GetName());
2531     if (branch) branch->SetAddress(&fTrackReferences);
2532   }
2533
2534 }
2535
2536 void AliTPC::FinishPrimary()
2537 {
2538   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();
2539   if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();
2540 }
2541
2542
2543 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2544 {
2545   //
2546   // add hit to the list
2547   Int_t rtrack;
2548   if (fIshunt) {
2549     int primary = gAlice->GetMCApp()->GetPrimary(track);
2550     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2551     rtrack=primary;
2552   } else {
2553     rtrack=track;
2554     gAlice->GetMCApp()->FlagTrack(track);
2555   }
2556   //AliTPChit *hit = (AliTPChit*)fHits->UncheckedAt(fNhits-1);
2557   //if (hit->fTrack!=rtrack)
2558   //  cout<<"bad track number\n";
2559   if (fTrackHits && fHitType&4)
2560     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2561                              hits[1],hits[2],(Int_t)hits[3]);
2562   if (fTrackHitsOld &&fHitType&2 )
2563     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2564                                 hits[1],hits[2],(Int_t)hits[3]);
2565
2566 }
2567
2568 void AliTPC::ResetHits()
2569 {
2570   if (fHitType&1) AliDetector::ResetHits();
2571   if (fHitType>1) ResetHits2();
2572 }
2573
2574 void AliTPC::ResetHits2()
2575 {
2576   //
2577   //reset hits
2578   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2579   if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2580
2581 }
2582
2583 AliHit* AliTPC::FirstHit(Int_t track)
2584 {
2585   if (fHitType>1) return FirstHit2(track);
2586   return AliDetector::FirstHit(track);
2587 }
2588 AliHit* AliTPC::NextHit()
2589 {
2590   //
2591   // gets next hit
2592   //
2593   if (fHitType>1) return NextHit2();
2594
2595   return AliDetector::NextHit();
2596 }
2597
2598 AliHit* AliTPC::FirstHit2(Int_t track)
2599 {
2600   //
2601   // Initialise the hit iterator
2602   // Return the address of the first hit for track
2603   // If track>=0 the track is read from disk
2604   // while if track<0 the first hit of the current
2605   // track is returned
2606   //
2607   if(track>=0) {
2608     gAlice->ResetHits();
2609     TreeH()->GetEvent(track);
2610   }
2611   //
2612   if (fTrackHits && fHitType&4) {
2613     fTrackHits->First();
2614     return fTrackHits->GetHit();
2615   }
2616   if (fTrackHitsOld && fHitType&2) {
2617     fTrackHitsOld->First();
2618     return fTrackHitsOld->GetHit();
2619   }
2620
2621   else return 0;
2622 }
2623
2624 AliHit* AliTPC::NextHit2()
2625 {
2626   //
2627   //Return the next hit for the current track
2628
2629
2630   if (fTrackHitsOld && fHitType&2) {
2631     fTrackHitsOld->Next();
2632     return fTrackHitsOld->GetHit();
2633   }
2634   if (fTrackHits) {
2635     fTrackHits->Next();
2636     return fTrackHits->GetHit();
2637   }
2638   else
2639     return 0;
2640 }
2641
2642 void AliTPC::LoadPoints(Int_t)
2643 {
2644   //
2645   Int_t a = 0;
2646   /*  if(fHitType==1) return AliDetector::LoadPoints(a);
2647   LoadPoints2(a);
2648   */
2649   if(fHitType==1) AliDetector::LoadPoints(a);
2650   else LoadPoints2(a);
2651
2652   // LoadPoints3(a);
2653
2654 }
2655
2656
2657 void AliTPC::RemapTrackHitIDs(Int_t *map)
2658 {
2659   //
2660   // remapping
2661   //
2662   if (!fTrackHits) return;
2663
2664   if (fTrackHitsOld && fHitType&2){
2665     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2666     for (UInt_t i=0;i<arr->GetSize();i++){
2667       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2668       info->fTrackID = map[info->fTrackID];
2669     }
2670   }
2671   if (fTrackHitsOld && fHitType&4){
2672     TClonesArray * arr = fTrackHits->GetArray();;
2673     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2674       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2675       info->fTrackID = map[info->fTrackID];
2676     }
2677   }
2678 }
2679
2680 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2681 {
2682   //return bool information - is track in given volume
2683   //load only part of the track information
2684   //return true if current track is in volume
2685   //
2686   //  return kTRUE;
2687   if (fTrackHitsOld && fHitType&2) {
2688     TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
2689     br->GetEvent(track);
2690     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2691     for (UInt_t j=0;j<ar->GetSize();j++){
2692       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2693     }
2694   }
2695
2696   if (fTrackHits && fHitType&4) {
2697     TBranch * br1 = TreeH()->GetBranch("fVolumes");
2698     TBranch * br2 = TreeH()->GetBranch("fNVolumes");
2699     br2->GetEvent(track);
2700     br1->GetEvent(track);
2701     Int_t *volumes = fTrackHits->GetVolumes();
2702     Int_t nvolumes = fTrackHits->GetNVolumes();
2703     if (!volumes && nvolumes>0) {
2704       printf("Problematic track\t%d\t%d",track,nvolumes);
2705       return kFALSE;
2706     }
2707     for (Int_t j=0;j<nvolumes; j++)
2708       if (volumes[j]==id) return kTRUE;
2709   }
2710
2711   if (fHitType&1) {
2712     TBranch * br = TreeH()->GetBranch("fSector");
2713     br->GetEvent(track);
2714     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2715       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2716     }
2717   }
2718   return kFALSE;
2719
2720 }
2721
2722 //_____________________________________________________________________________
2723 void AliTPC::LoadPoints2(Int_t)
2724 {
2725   //
2726   // Store x, y, z of all hits in memory
2727   //
2728   if (fTrackHits == 0 && fTrackHitsOld==0) return;
2729   //
2730   Int_t nhits =0;
2731   if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
2732   if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
2733
2734   if (nhits == 0) return;
2735   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2736   if (fPoints == 0) fPoints = new TObjArray(tracks);
2737   AliHit *ahit;
2738   //
2739   Int_t *ntrk=new Int_t[tracks];
2740   Int_t *limi=new Int_t[tracks];
2741   Float_t **coor=new Float_t*[tracks];
2742   for(Int_t i=0;i<tracks;i++) {
2743     ntrk[i]=0;
2744     coor[i]=0;
2745     limi[i]=0;
2746   }
2747   //
2748   AliPoints *points = 0;
2749   Float_t *fp=0;
2750   Int_t trk;
2751   Int_t chunk=nhits/4+1;
2752   //
2753   // Loop over all the hits and store their position
2754   //
2755   ahit = FirstHit2(-1);
2756   while (ahit){
2757     trk=ahit->GetTrack();
2758     if(ntrk[trk]==limi[trk]) {
2759       //
2760       // Initialise a new track
2761       fp=new Float_t[3*(limi[trk]+chunk)];
2762       if(coor[trk]) {
2763         memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2764         delete [] coor[trk];
2765       }
2766       limi[trk]+=chunk;
2767       coor[trk] = fp;
2768     } else {
2769       fp = coor[trk];
2770     }
2771     fp[3*ntrk[trk]  ] = ahit->X();
2772     fp[3*ntrk[trk]+1] = ahit->Y();
2773     fp[3*ntrk[trk]+2] = ahit->Z();
2774     ntrk[trk]++;
2775     ahit = NextHit2();
2776   }
2777
2778
2779
2780   //
2781   for(trk=0; trk<tracks; ++trk) {
2782     if(ntrk[trk]) {
2783       points = new AliPoints();
2784       points->SetMarkerColor(GetMarkerColor());
2785       points->SetMarkerSize(GetMarkerSize());
2786       points->SetDetector(this);
2787       points->SetParticle(trk);
2788       points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
2789       fPoints->AddAt(points,trk);
2790       delete [] coor[trk];
2791       coor[trk]=0;
2792     }
2793   }
2794   delete [] coor;
2795   delete [] ntrk;
2796   delete [] limi;
2797 }
2798
2799
2800 //_____________________________________________________________________________
2801 void AliTPC::LoadPoints3(Int_t)
2802 {
2803   //
2804   // Store x, y, z of all hits in memory
2805   // - only intersection point with pad row
2806   if (fTrackHits == 0) return;
2807   //
2808   Int_t nhits = fTrackHits->GetEntriesFast();
2809   if (nhits == 0) return;
2810   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2811   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2812   fPoints->Expand(2*tracks);
2813   AliHit *ahit;
2814   //
2815   Int_t *ntrk=new Int_t[tracks];
2816   Int_t *limi=new Int_t[tracks];
2817   Float_t **coor=new Float_t*[tracks];
2818   for(Int_t i=0;i<tracks;i++) {
2819     ntrk[i]=0;
2820     coor[i]=0;
2821     limi[i]=0;
2822   }
2823   //
2824   AliPoints *points = 0;
2825   Float_t *fp=0;
2826   Int_t trk;
2827   Int_t chunk=nhits/4+1;
2828   //
2829   // Loop over all the hits and store their position
2830   //
2831   ahit = FirstHit2(-1);
2832   //for (Int_t hit=0;hit<nhits;hit++) {
2833
2834   Int_t lastrow = -1;
2835   while (ahit){
2836     //    ahit = (AliHit*)fHits->UncheckedAt(hit);
2837     trk=ahit->GetTrack();
2838     Float_t  x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2839     Int_t    index[3]={1,((AliTPChit*)ahit)->fSector,0};
2840     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2841     if (currentrow!=lastrow){
2842       lastrow = currentrow;
2843       //later calculate intersection point
2844       if(ntrk[trk]==limi[trk]) {
2845         //
2846         // Initialise a new track
2847         fp=new Float_t[3*(limi[trk]+chunk)];
2848         if(coor[trk]) {
2849           memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2850           delete [] coor[trk];
2851         }
2852         limi[trk]+=chunk;
2853         coor[trk] = fp;
2854       } else {
2855         fp = coor[trk];
2856       }
2857       fp[3*ntrk[trk]  ] = ahit->X();
2858       fp[3*ntrk[trk]+1] = ahit->Y();
2859       fp[3*ntrk[trk]+2] = ahit->Z();
2860       ntrk[trk]++;
2861     }
2862     ahit = NextHit2();
2863   }
2864
2865   //
2866   for(trk=0; trk<tracks; ++trk) {
2867     if(ntrk[trk]) {
2868       points = new AliPoints();
2869       points->SetMarkerColor(GetMarkerColor()+1);
2870       points->SetMarkerStyle(5);
2871       points->SetMarkerSize(0.2);
2872       points->SetDetector(this);
2873       points->SetParticle(trk);
2874       //      points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle()20);
2875       points->SetPolyMarker(ntrk[trk],coor[trk],30);
2876       fPoints->AddAt(points,tracks+trk);
2877       delete [] coor[trk];
2878       coor[trk]=0;
2879     }
2880   }
2881   delete [] coor;
2882   delete [] ntrk;
2883   delete [] limi;
2884 }
2885
2886
2887
2888 void AliTPC::FindTrackHitsIntersection(TClonesArray * /*arr*/)
2889 {
2890
2891   //
2892   //fill clones array with intersection of current point with the
2893   //middle of the row
2894   Int_t sector;
2895   Int_t ipart;
2896
2897   const Int_t kcmaxhits=30000;
2898   AliTPCFastVector * xxxx = new AliTPCFastVector(kcmaxhits*4);
2899   AliTPCFastVector & xxx = *xxxx;
2900   Int_t maxhits = kcmaxhits;
2901
2902   //
2903   AliTPChit * tpcHit=0;
2904   tpcHit = (AliTPChit*)FirstHit2(-1);
2905   Int_t currentIndex=0;
2906   Int_t lastrow=-1;  //last writen row
2907
2908   while (tpcHit){
2909     if (tpcHit==0) continue;
2910     sector=tpcHit->fSector; // sector number
2911     ipart=tpcHit->Track();
2912
2913     //find row number
2914
2915     Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
2916     Int_t    index[3]={1,sector,0};
2917     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2918     if (currentrow<0) continue;
2919     if (lastrow<0) lastrow=currentrow;
2920     if (currentrow==lastrow){
2921       if ( currentIndex>=maxhits){
2922         maxhits+=kcmaxhits;
2923         xxx.ResizeTo(4*maxhits);
2924       }
2925       xxx(currentIndex*4)=x[0];
2926       xxx(currentIndex*4+1)=x[1];
2927       xxx(currentIndex*4+2)=x[2];
2928       xxx(currentIndex*4+3)=tpcHit->fQ;
2929       currentIndex++;
2930     }
2931     else
2932       if (currentIndex>2){
2933         Float_t sumx=0;
2934         Float_t sumx2=0;
2935         Float_t sumx3=0;
2936         Float_t sumx4=0;
2937         Float_t sumy=0;
2938         Float_t sumxy=0;
2939         Float_t sumx2y=0;
2940         Float_t sumz=0;
2941         Float_t sumxz=0;
2942         Float_t sumx2z=0;
2943         Float_t sumq=0;
2944         for (Int_t index=0;index<currentIndex;index++){
2945           Float_t x,x2,x3,x4;
2946           x=x2=x3=x4=xxx(index*4);
2947           x2*=x;
2948           x3*=x2;
2949           x4*=x3;
2950           sumx+=x;
2951           sumx2+=x2;
2952           sumx3+=x3;
2953           sumx4+=x4;
2954           sumy+=xxx(index*4+1);
2955           sumxy+=xxx(index*4+1)*x;
2956           sumx2y+=xxx(index*4+1)*x2;
2957           sumz+=xxx(index*4+2);
2958           sumxz+=xxx(index*4+2)*x;
2959           sumx2z+=xxx(index*4+2)*x2;
2960           sumq+=xxx(index*4+3);
2961         }
2962         Float_t centralPad = (fTPCParam->GetNPads(sector,lastrow)-1)/2;
2963         Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
2964           sumx2*(sumx*sumx3-sumx2*sumx2);
2965
2966         Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
2967           sumx2*(sumxy*sumx3-sumx2y*sumx2);
2968         Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
2969           sumx2*(sumxz*sumx3-sumx2z*sumx2);
2970
2971         Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
2972           sumx2*(sumx*sumx2y-sumx2*sumxy);
2973         Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
2974           sumx2*(sumx*sumx2z-sumx2*sumxz);
2975
2976         Float_t y=detay/det+centralPad;
2977         Float_t z=detaz/det;
2978         Float_t by=detby/det; //y angle
2979         Float_t bz=detbz/det; //z angle
2980         sumy/=Float_t(currentIndex);
2981         sumz/=Float_t(currentIndex);
2982
2983         AliComplexCluster cl;
2984         cl.fX=z;
2985         cl.fY=y;
2986         cl.fQ=sumq;
2987         cl.fSigmaX2=bz;
2988         cl.fSigmaY2=by;
2989         cl.fTracks[0]=ipart;
2990
2991         AliTPCClustersRow * row = (fClustersArray->GetRow(sector,lastrow));
2992         if (row!=0) row->InsertCluster(&cl);
2993         currentIndex=0;
2994         lastrow=currentrow;
2995       } //end of calculating cluster for given row
2996
2997   } // end of loop over hits
2998   xxxx->Delete();
2999
3000 }
3001 //_______________________________________________________________________________
3002 void AliTPC::Digits2Reco(Int_t firstevent,Int_t lastevent)
3003 {
3004   // produces rec points from digits and writes them on the root file
3005   // AliTPCclusters.root
3006
3007   TDirectory *cwd = gDirectory;
3008
3009
3010   AliTPCParamSR *dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60");
3011   if(dig){
3012     printf("You are running 2 pad-length geom hits with 3 pad-length geom digits\n");
3013     delete dig;
3014     dig = new AliTPCParamSR();
3015   }
3016   else
3017   {
3018    dig=(AliTPCParamSR *)gDirectory->Get("75x40_100x60_150x60");
3019   }
3020   if(!dig){
3021    printf("No TPC parameters found\n");
3022    exit(3);
3023   }
3024
3025   SetParam(dig);
3026   cout<<"AliTPC::Digits2Reco: TPC parameteres have been set"<<endl;
3027   TFile *out;
3028   if(!gSystem->Getenv("CONFIG_FILE")){
3029     out=TFile::Open("AliTPCclusters.root","recreate");
3030   }
3031   else {
3032     const char *tmp1;
3033     const char *tmp2;
3034     char tmp3[80];
3035     tmp1=gSystem->Getenv("CONFIG_FILE_PREFIX");
3036     tmp2=gSystem->Getenv("CONFIG_OUTDIR");
3037     sprintf(tmp3,"%s%s/AliTPCclusters.root",tmp1,tmp2);
3038     out=TFile::Open(tmp3,"recreate");
3039   }
3040
3041   TStopwatch timer;
3042   cout<<"AliTPC::Digits2Reco - determination of rec points begins"<<endl;
3043   timer.Start();
3044   cwd->cd();
3045   for(Int_t iev=firstevent;iev<lastevent+1;iev++){
3046
3047     printf("Processing event %d\n",iev);
3048     Digits2Clusters(iev);
3049   }
3050   cout<<"AliTPC::Digits2Reco - determination of rec points ended"<<endl;
3051   timer.Stop();
3052   timer.Print();
3053   out->Close();
3054   cwd->cd();
3055
3056
3057 }
3058
3059 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
3060 {
3061 //Makes TPC loader
3062  fLoader = new AliTPCLoader(GetName(),topfoldername);
3063  return fLoader;
3064 }
3065
3066 ////////////////////////////////////////////////////////////////////////
3067 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
3068 //
3069 // load TPC paarmeters from a given file or create new if the object
3070 // is not found there
3071 // 12/05/2003 This method should be moved to the AliTPCLoader
3072 // and one has to decide where to store the TPC parameters
3073 // M.Kowalski
3074   char paramName[50];
3075   sprintf(paramName,"75x40_100x60_150x60");
3076   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
3077   if (paramTPC) {
3078     cout<<"TPC parameters "<<paramName<<" found."<<endl;
3079   } else {
3080     cerr<<"TPC parameters not found. Create new (they may be incorrect)."
3081         <<endl;
3082     paramTPC = new AliTPCParamSR;
3083   }
3084   return paramTPC;
3085
3086 // the older version of parameters can be accessed with this code.
3087 // In some cases, we have old parameters saved in the file but
3088 // digits were created with new parameters, it can be distinguish
3089 // by the name of TPC TreeD. The code here is just for the case
3090 // we would need to compare with old data, uncomment it if needed.
3091 //
3092 //  char paramName[50];
3093 //  sprintf(paramName,"75x40_100x60");
3094 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
3095 //  if (paramTPC) {
3096 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
3097 //  } else {
3098 //    sprintf(paramName,"75x40_100x60_150x60");
3099 //    paramTPC=(AliTPCParam*)in->Get(paramName);
3100 //    if (paramTPC) {
3101 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
3102 //    } else {
3103 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
3104 //          <<endl;
3105 //      paramTPC = new AliTPCParamSR;
3106 //    }
3107 //  }
3108 //  return paramTPC;
3109
3110 }
3111
3112
3113 //____________________________________________________________________________
3114 Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
3115 {
3116   //
3117   // Parametrised error of the cluster reconstruction (pad direction)
3118   //
3119   // Sigma rphi
3120   const Float_t kArphi=0.41818e-2;
3121   const Float_t kBrphi=0.17460e-4;
3122   const Float_t kCrphi=0.30993e-2;
3123   const Float_t kDrphi=0.41061e-3;
3124
3125   pt=TMath::Abs(pt)*1000.;
3126   Double_t x=r/pt;
3127   tgl=TMath::Abs(tgl);
3128   Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
3129   if (s<0.4e-3) s=0.4e-3;
3130   s*=1.3; //Iouri Belikov
3131
3132   return s;
3133 }
3134
3135
3136 //____________________________________________________________________________
3137 Double_t SigmaZ2(Double_t r, Double_t tgl)
3138 {
3139   //
3140   // Parametrised error of the cluster reconstruction (drift direction)
3141   //
3142   // Sigma z
3143   const Float_t kAz=0.39614e-2;
3144   const Float_t kBz=0.22443e-4;
3145   const Float_t kCz=0.51504e-1;
3146
3147
3148   tgl=TMath::Abs(tgl);
3149   Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
3150   if (s<0.4e-3) s=0.4e-3;
3151   s*=1.3; //Iouri Belikov
3152
3153   return s;
3154 }
3155