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 <TVirtualMC.h>
  53 #include <TString.h>
  54 #include <TF2.h>
  55 #include <TStopwatch.h>
  56
  57 #include "AliArrayBranch.h"
  58 #include "AliDigits.h"
  59 #include "AliMagF.h"
  60 #include "AliPoints.h"
  61 #include "AliRun.h"
  62 #include "AliRunLoader.h"
  63 #include "AliSimDigits.h"
  64 #include "AliTPC.h"
  65 #include "AliTPC.h"
  66 #include "AliTPCDigitsArray.h"
  67 #include "AliTPCLoader.h"
  68 #include "AliTPCPRF2D.h"
  69 #include "AliTPCParamSR.h"
  70 #include "AliTPCRF1D.h"
  71 #include "AliTPCTrackHits.h"
  72 #include "AliTPCTrackHitsV2.h"
  73 #include "AliTrackReference.h"
  74 #include "AliMC.h"
  75 #include "AliTPCDigitizer.h"
  76 #include "AliTPCBuffer.h"
  77 #include "AliTPCDDLRawData.h"
  78 #include "AliLog.h"
  79
  80
  81 ClassImp(AliTPC)
  82 //_____________________________________________________________________________
  83 AliTPC::AliTPC()
  84 {
  85   //
  86   // Default constructor
  87   //
  88   fIshunt   = 0;
  89   fHits     = 0;
  90   fDigits   = 0;
  91   fNsectors = 0;
  92   fDigitsArray = 0;
  93   fDefaults = 0;
  94   fTrackHits = 0;
  95   fTrackHitsOld = 0;
  96 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
  97   fHitType = 4; // ROOT containers
  98 #else
  99   fHitType = 2; //default CONTAINERS - based on ROOT structure
 100 #endif
 101   fTPCParam = 0;
 102   fNoiseTable = 0;
 103   fActiveSectors =0;
 104
 105 }
 106
 107 //_____________________________________________________________________________
 108 AliTPC::AliTPC(const char *name, const char *title)
 109       : AliDetector(name,title)
 110 {
 111   //
 112   // Standard constructor
 113   //
 114
 115   //
 116   // Initialise arrays of hits and digits
 117   fHits     = new TClonesArray("AliTPChit",  176);
 118   gAlice->GetMCApp()->AddHitList(fHits);
 119   fDigitsArray = 0;
 120   fDefaults = 0;
 121   //
 122   fTrackHits = new AliTPCTrackHitsV2;
 123   fTrackHits->SetHitPrecision(0.002);
 124   fTrackHits->SetStepPrecision(0.003);
 125   fTrackHits->SetMaxDistance(100);
 126
 127   fTrackHitsOld = new AliTPCTrackHits;  //MI - 13.09.2000
 128   fTrackHitsOld->SetHitPrecision(0.002);
 129   fTrackHitsOld->SetStepPrecision(0.003);
 130   fTrackHitsOld->SetMaxDistance(100);
 131
 132   fNoiseTable =0;
 133
 134 #if ROOT_VERSION_CODE >= ROOT_VERSION(4,0,1)
 135   fHitType = 4; // ROOT containers
 136 #else
 137   fHitType = 2;
 138 #endif
 139   fActiveSectors = 0;
 140   //
 141   // Initialise counters
 142   fNsectors = 0;
 143
 144   //
 145   fIshunt     =  0;
 146   //
 147   // Initialise color attributes
 148   SetMarkerColor(kYellow);
 149
 150   //
 151   //  Set TPC parameters
 152   //
 153
 154
 155   if (!strcmp(title,"Default")) {
 156     fTPCParam = new AliTPCParamSR;
 157   } else {
 158     AliWarning("In Config.C you must set non-default parameters.");
 159     fTPCParam=0;
 160   }
 161
 162 }
 163
 164 //_____________________________________________________________________________
 165 AliTPC::AliTPC(const AliTPC& t):AliDetector(t){
 166   //
 167   // dummy copy constructor
 168   //
 169 }
 170 AliTPC::~AliTPC()
 171 {
 172   //
 173   // TPC destructor
 174   //
 175
 176   fIshunt   = 0;
 177   delete fHits;
 178   delete fDigits;
 179   delete fTPCParam;
 180   delete fTrackHits; //MI 15.09.2000
 181   delete fTrackHitsOld; //MI 10.12.2001
 182   if (fNoiseTable) delete [] fNoiseTable;
 183
 184 }
 185
 186 //_____________________________________________________________________________
 187 void AliTPC::AddHit(Int_t track, Int_t *vol, Float_t *hits)
 188 {
 189   //
 190   // Add a hit to the list
 191   //
 192   if (fHitType&1){
 193     TClonesArray &lhits = *fHits;
 194     new(lhits[fNhits++]) AliTPChit(fIshunt,track,vol,hits);
 195   }
 196   if (fHitType>1)
 197     AddHit2(track,vol,hits);
 198 }
 199
 200 //_____________________________________________________________________________
 201 void AliTPC::BuildGeometry()
 202 {
 203
 204   //
 205   // Build TPC ROOT TNode geometry for the event display
 206   //
 207   TNode *nNode, *nTop;
 208   TTUBS *tubs;
 209   Int_t i;
 210   const int kColorTPC=19;
 211   char name[5], title[25];
 212   const Double_t kDegrad=TMath::Pi()/180;
 213   const Double_t kRaddeg=180./TMath::Pi();
 214
 215
 216   Float_t innerOpenAngle = fTPCParam->GetInnerAngle();
 217   Float_t outerOpenAngle = fTPCParam->GetOuterAngle();
 218
 219   Float_t innerAngleShift = fTPCParam->GetInnerAngleShift();
 220   Float_t outerAngleShift = fTPCParam->GetOuterAngleShift();
 221
 222   Int_t nLo = fTPCParam->GetNInnerSector()/2;
 223   Int_t nHi = fTPCParam->GetNOuterSector()/2;
 224
 225   const Double_t kloAng = (Double_t)TMath::Nint(innerOpenAngle*kRaddeg);
 226   const Double_t khiAng = (Double_t)TMath::Nint(outerOpenAngle*kRaddeg);
 227   const Double_t kloAngSh = (Double_t)TMath::Nint(innerAngleShift*kRaddeg);
 228   const Double_t khiAngSh = (Double_t)TMath::Nint(outerAngleShift*kRaddeg);
 229
 230
 231   const Double_t kloCorr = 1/TMath::Cos(0.5*kloAng*kDegrad);
 232   const Double_t khiCorr = 1/TMath::Cos(0.5*khiAng*kDegrad);
 233
 234   Double_t rl,ru;
 235
 236
 237   //
 238   // Get ALICE top node
 239   //
 240
 241   nTop=gAlice->GetGeometry()->GetNode("alice");
 242
 243   //  inner sectors
 244
 245   rl = fTPCParam->GetInnerRadiusLow();
 246   ru = fTPCParam->GetInnerRadiusUp();
 247
 248
 249   for(i=0;i<nLo;i++) {
 250     sprintf(name,"LS%2.2d",i);
 251     name[4]='\0';
 252     sprintf(title,"TPC low sector %3d",i);
 253     title[24]='\0';
 254
 255     tubs = new TTUBS(name,title,"void",rl*kloCorr,ru*kloCorr,250.,
 256                      kloAng*(i-0.5)+kloAngSh,kloAng*(i+0.5)+kloAngSh);
 257     tubs->SetNumberOfDivisions(1);
 258     nTop->cd();
 259     nNode = new TNode(name,title,name,0,0,0,"");
 260     nNode->SetLineColor(kColorTPC);
 261     fNodes->Add(nNode);
 262   }
 263
 264   // Outer sectors
 265
 266   rl = fTPCParam->GetOuterRadiusLow();
 267   ru = fTPCParam->GetOuterRadiusUp();
 268
 269   for(i=0;i<nHi;i++) {
 270     sprintf(name,"US%2.2d",i);
 271     name[4]='\0';
 272     sprintf(title,"TPC upper sector %d",i);
 273     title[24]='\0';
 274     tubs = new TTUBS(name,title,"void",rl*khiCorr,ru*khiCorr,250,
 275                      khiAng*(i-0.5)+khiAngSh,khiAng*(i+0.5)+khiAngSh);
 276     tubs->SetNumberOfDivisions(1);
 277     nTop->cd();
 278     nNode = new TNode(name,title,name,0,0,0,"");
 279     nNode->SetLineColor(kColorTPC);
 280     fNodes->Add(nNode);
 281   }
 282
 283 }
 284
 285 //_____________________________________________________________________________
 286 Int_t AliTPC::DistancetoPrimitive(Int_t , Int_t ) const
 287 {
 288   //
 289   // Calculate distance from TPC to mouse on the display
 290   // Dummy procedure
 291   //
 292   return 9999;
 293 }
 294
 295
 296 //_____________________________________________________________________________
 297 void AliTPC::CreateMaterials()
 298 {
 299   //-----------------------------------------------
 300   // Create Materials for for TPC simulations
 301   //-----------------------------------------------
 302
 303   //-----------------------------------------------------------------
 304   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
 305   //-----------------------------------------------------------------
 306
 307   Int_t iSXFLD=gAlice->Field()->Integ();
 308   Float_t sXMGMX=gAlice->Field()->Max();
 309
 310   Float_t amat[5]; // atomic numbers
 311   Float_t zmat[5]; // z
 312   Float_t wmat[5]; // proportions
 313
 314   Float_t density;
 315   Float_t apure[2];
 316
 317
 318   //***************** Gases *************************
 319
 320   //-------------------------------------------------
 321   // pure gases
 322   //-------------------------------------------------
 323
 324   // Neon
 325
 326
 327   amat[0]= 20.18;
 328   zmat[0]= 10.;
 329   density = 0.0009;
 330
 331   apure[0]=amat[0];
 332
 333   AliMaterial(20,"Ne",amat[0],zmat[0],density,999.,999.);
 334
 335   // Argon
 336
 337   amat[0]= 39.948;
 338   zmat[0]= 18.;
 339   density = 0.001782;
 340
 341   apure[1]=amat[0];
 342
 343   AliMaterial(21,"Ar",amat[0],zmat[0],density,999.,999.);
 344
 345
 346   //--------------------------------------------------------------
 347   // gases - compounds
 348   //--------------------------------------------------------------
 349
 350   Float_t amol[3];
 351
 352   // CO2
 353
 354   amat[0]=12.011;
 355   amat[1]=15.9994;
 356
 357   zmat[0]=6.;
 358   zmat[1]=8.;
 359
 360   wmat[0]=1.;
 361   wmat[1]=2.;
 362
 363   density=0.001977;
 364
 365   amol[0] = amat[0]*wmat[0]+amat[1]*wmat[1];
 366
 367   AliMixture(10,"CO2",amat,zmat,density,-2,wmat);
 368
 369   // CF4
 370
 371   amat[0]=12.011;
 372   amat[1]=18.998;
 373
 374   zmat[0]=6.;
 375   zmat[1]=9.;
 376
 377   wmat[0]=1.;
 378   wmat[1]=4.;
 379
 380   density=0.003034;
 381
 382   amol[1] = amat[0]*wmat[0]+amat[1]*wmat[1];
 383
 384   AliMixture(11,"CF4",amat,zmat,density,-2,wmat);
 385
 386
 387   // CH4
 388
 389   amat[0]=12.011;
 390   amat[1]=1.;
 391
 392   zmat[0]=6.;
 393   zmat[1]=1.;
 394
 395   wmat[0]=1.;
 396   wmat[1]=4.;
 397
 398   density=0.000717;
 399
 400   amol[2] = amat[0]*wmat[0]+amat[1]*wmat[1];
 401
 402   AliMixture(12,"CH4",amat,zmat,density,-2,wmat);
 403
 404   //----------------------------------------------------------------
 405   // gases - mixtures, ID >= 20 pure gases, <= 10 ID < 20 -compounds
 406   //----------------------------------------------------------------
 407
 408   char namate[21]="";
 409   density = 0.;
 410   Float_t am=0;
 411   Int_t nc;
 412   Float_t rho,absl,x0,buf[1];
 413   Int_t nbuf;
 414   Float_t a,z;
 415
 416   for(nc = 0;nc<fNoComp;nc++) {
 417
 418     // retrive material constants
 419
 420     gMC->Gfmate((*fIdmate)[fMixtComp[nc]],namate,a,z,rho,x0,absl,buf,nbuf);
 421
 422     amat[nc] = a;
 423     zmat[nc] = z;
 424
 425     Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 426
 427     am += fMixtProp[nc]*((fMixtComp[nc]>=20) ? apure[nnc] : amol[nnc]);
 428     density += fMixtProp[nc]*rho;  // density of the mixture
 429
 430   }
 431
 432   // mixture proportions by weight!
 433
 434   for(nc = 0;nc<fNoComp;nc++) {
 435
 436     Int_t nnc = (fMixtComp[nc]>=20) ? fMixtComp[nc]%20 : fMixtComp[nc]%10;
 437
 438     wmat[nc] = fMixtProp[nc]*((fMixtComp[nc]>=20) ?
 439                               apure[nnc] : amol[nnc])/am;
 440
 441   }
 442
 443   // Drift gases 1 - nonsensitive, 2 - sensitive
 444
 445   AliMixture(31,"Drift gas 1",amat,zmat,density,fNoComp,wmat);
 446   AliMixture(32,"Drift gas 2",amat,zmat,density,fNoComp,wmat);
 447
 448
 449   // Air
 450
 451   amat[0] = 14.61;
 452   zmat[0] = 7.3;
 453   density = 0.001205;
 454
 455   AliMaterial(24,"Air",amat[0],zmat[0],density,999.,999.);
 456
 457
 458   //----------------------------------------------------------------------
 459   //               solid materials
 460   //----------------------------------------------------------------------
 461
 462
 463   // Kevlar C14H22O2N2
 464
 465   amat[0] = 12.011;
 466   amat[1] = 1.;
 467   amat[2] = 15.999;
 468   amat[3] = 14.006;
 469
 470   zmat[0] = 6.;
 471   zmat[1] = 1.;
 472   zmat[2] = 8.;
 473   zmat[3] = 7.;
 474
 475   wmat[0] = 14.;
 476   wmat[1] = 22.;
 477   wmat[2] = 2.;
 478   wmat[3] = 2.;
 479
 480   density = 1.45;
 481
 482   AliMixture(34,"Kevlar",amat,zmat,density,-4,wmat);
 483
 484   // NOMEX
 485
 486   amat[0] = 12.011;
 487   amat[1] = 1.;
 488   amat[2] = 15.999;
 489   amat[3] = 14.006;
 490
 491   zmat[0] = 6.;
 492   zmat[1] = 1.;
 493   zmat[2] = 8.;
 494   zmat[3] = 7.;
 495
 496   wmat[0] = 14.;
 497   wmat[1] = 22.;
 498   wmat[2] = 2.;
 499   wmat[3] = 2.;
 500
 501   density = 0.03;
 502
 503
 504   AliMixture(35,"NOMEX",amat,zmat,density,-4,wmat);
 505
 506   // Makrolon C16H18O3
 507
 508   amat[0] = 12.011;
 509   amat[1] = 1.;
 510   amat[2] = 15.999;
 511
 512   zmat[0] = 6.;
 513   zmat[1] = 1.;
 514   zmat[2] = 8.;
 515
 516   wmat[0] = 16.;
 517   wmat[1] = 18.;
 518   wmat[2] = 3.;
 519
 520   density = 1.2;
 521
 522   AliMixture(36,"Makrolon",amat,zmat,density,-3,wmat);
 523
 524   // Mylar C5H4O2
 525
 526   amat[0]=12.011;
 527   amat[1]=1.;
 528   amat[2]=15.9994;
 529
 530   zmat[0]=6.;
 531   zmat[1]=1.;
 532   zmat[2]=8.;
 533
 534   wmat[0]=5.;
 535   wmat[1]=4.;
 536   wmat[2]=2.;
 537
 538   density = 1.39;
 539
 540   AliMixture(37, "Mylar",amat,zmat,density,-3,wmat);
 541
 542   // SiO2 - used later for the glass fiber
 543
 544   amat[0]=28.086;
 545   amat[1]=15.9994;
 546
 547   zmat[0]=14.;
 548   zmat[1]=8.;
 549
 550   wmat[0]=1.;
 551   wmat[1]=2.;
 552
 553
 554   AliMixture(38,"SiO2",amat,zmat,2.2,-2,wmat); //SiO2 - quartz (rho=2.2)
 555
 556   // Al
 557
 558   amat[0] = 26.98;
 559   zmat[0] = 13.;
 560
 561   density = 2.7;
 562
 563   AliMaterial(40,"Al",amat[0],zmat[0],density,999.,999.);
 564
 565   // Si
 566
 567   amat[0] = 28.086;
 568   zmat[0] = 14.;
 569
 570   density = 2.33;
 571
 572   AliMaterial(41,"Si",amat[0],zmat[0],density,999.,999.);
 573
 574   // Cu
 575
 576   amat[0] = 63.546;
 577   zmat[0] = 29.;
 578
 579   density = 8.96;
 580
 581   AliMaterial(42,"Cu",amat[0],zmat[0],density,999.,999.);
 582
 583   // Tedlar C2H3F
 584
 585   amat[0] = 12.011;
 586   amat[1] = 1.;
 587   amat[2] = 18.998;
 588
 589   zmat[0] = 6.;
 590   zmat[1] = 1.;
 591   zmat[2] = 9.;
 592
 593   wmat[0] = 2.;
 594   wmat[1] = 3.;
 595   wmat[2] = 1.;
 596
 597   density = 1.71;
 598
 599   AliMixture(43, "Tedlar",amat,zmat,density,-3,wmat);
 600
 601
 602   // Plexiglas  C5H8O2
 603
 604   amat[0]=12.011;
 605   amat[1]=1.;
 606   amat[2]=15.9994;
 607
 608   zmat[0]=6.;
 609   zmat[1]=1.;
 610   zmat[2]=8.;
 611
 612   wmat[0]=5.;
 613   wmat[1]=8.;
 614   wmat[2]=2.;
 615
 616   density=1.18;
 617
 618   AliMixture(44,"Plexiglas",amat,zmat,density,-3,wmat);
 619
 620   // Epoxy - C14 H20 O3
 621
 622
 623   amat[0]=12.011;
 624   amat[1]=1.;
 625   amat[2]=15.9994;
 626
 627   zmat[0]=6.;
 628   zmat[1]=1.;
 629   zmat[2]=8.;
 630
 631   wmat[0]=14.;
 632   wmat[1]=20.;
 633   wmat[2]=3.;
 634
 635   density=1.25;
 636
 637   AliMixture(45,"Epoxy",amat,zmat,density,-3,wmat);
 638
 639   // Carbon
 640
 641   amat[0]=12.011;
 642   zmat[0]=6.;
 643   density= 2.265;
 644
 645   AliMaterial(46,"C",amat[0],zmat[0],density,999.,999.);
 646
 647   // get epoxy
 648
 649   gMC->Gfmate((*fIdmate)[45],namate,amat[1],zmat[1],rho,x0,absl,buf,nbuf);
 650
 651   // Carbon fiber
 652
 653   wmat[0]=0.644; // by weight!
 654   wmat[1]=0.356;
 655
 656   density=0.5*(1.25+2.265);
 657
 658   AliMixture(47,"Cfiber",amat,zmat,density,2,wmat);
 659
 660   // get SiO2
 661
 662   gMC->Gfmate((*fIdmate)[38],namate,amat[0],zmat[0],rho,x0,absl,buf,nbuf);
 663
 664   wmat[0]=0.725; // by weight!
 665   wmat[1]=0.275;
 666
 667   density=1.7;
 668
 669   AliMixture(39,"G10",amat,zmat,density,2,wmat);
 670
 671
 672
 673
 674   //----------------------------------------------------------
 675   // tracking media for gases
 676   //----------------------------------------------------------
 677
 678   AliMedium(0, "Air", 24, 0, iSXFLD, sXMGMX, 10., 999., .1, .01, .1);
 679   AliMedium(1, "Drift gas 1", 31, 0, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 680   AliMedium(2, "Drift gas 2", 32, 1, iSXFLD, sXMGMX, 10., 999.,.1,.001, .001);
 681   AliMedium(3,"CO2",10,0, iSXFLD, sXMGMX, 10., 999.,.1, .001, .001);
 682
 683   //-----------------------------------------------------------
 684   // tracking media for solids
 685   //-----------------------------------------------------------
 686
 687   AliMedium(4,"Al",40,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 688   AliMedium(5,"Kevlar",34,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 689   AliMedium(6,"Nomex",35,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 690   AliMedium(7,"Makrolon",36,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 691   AliMedium(8,"Mylar",37,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 692   AliMedium(9,"Tedlar",43,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 693   AliMedium(10,"Cu",42,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 694   AliMedium(11,"Si",41,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 695   AliMedium(12,"G10",39,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 696   AliMedium(13,"Plexiglas",44,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 697   AliMedium(14,"Epoxy",45,0, iSXFLD, sXMGMX, 10., 999., .1, .0005, .001);
 698   AliMedium(15,"Cfiber",47,0, iSXFLD, sXMGMX, 10., 999., .1, .001, .001);
 699
 700 }
 701
 702 void AliTPC::GenerNoise(Int_t tablesize)
 703 {
 704   //
 705   //Generate table with noise
 706   //
 707   if (fTPCParam==0) {
 708     // error message
 709     fNoiseDepth=0;
 710     return;
 711   }
 712   if (fNoiseTable)  delete[] fNoiseTable;
 713   fNoiseTable = new Float_t[tablesize];
 714   fNoiseDepth = tablesize;
 715   fCurrentNoise =0; //!index of the noise in  the noise table
 716
 717   Float_t norm = fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac();
 718   for (Int_t i=0;i<tablesize;i++) fNoiseTable[i]= gRandom->Gaus(0,norm);
 719 }
 720
 721 Float_t AliTPC::GetNoise()
 722 {
 723   // get noise from table
 724   //  if ((fCurrentNoise%10)==0)
 725   //  fCurrentNoise= gRandom->Rndm()*fNoiseDepth;
 726   if (fCurrentNoise>=fNoiseDepth) fCurrentNoise=0;
 727   return fNoiseTable[fCurrentNoise++];
 728   //gRandom->Gaus(0, fTPCParam->GetNoise()*fTPCParam->GetNoiseNormFac());
 729 }
 730
 731
 732 Bool_t  AliTPC::IsSectorActive(Int_t sec) const
 733 {
 734   //
 735   // check if the sector is active
 736   if (!fActiveSectors) return kTRUE;
 737   else return fActiveSectors[sec];
 738 }
 739
 740 void    AliTPC::SetActiveSectors(Int_t * sectors, Int_t n)
 741 {
 742   // activate interesting sectors
 743   SetTreeAddress();//just for security
 744   if (fActiveSectors) delete [] fActiveSectors;
 745   fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 746   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 747   for (Int_t i=0;i<n;i++)
 748     if ((sectors[i]>=0) && sectors[i]<fTPCParam->GetNSector())  fActiveSectors[sectors[i]]=kTRUE;
 749
 750 }
 751
 752 void    AliTPC::SetActiveSectors(Int_t flag)
 753 {
 754   //
 755   // activate sectors which were hitted by tracks
 756   //loop over tracks
 757   SetTreeAddress();//just for security
 758   if (fHitType==0) return;  // if Clones hit - not short volume ID information
 759   if (fActiveSectors) delete [] fActiveSectors;
 760   fActiveSectors = new Bool_t[fTPCParam->GetNSector()];
 761   if (flag) {
 762     for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kTRUE;
 763     return;
 764   }
 765   for (Int_t i=0;i<fTPCParam->GetNSector();i++) fActiveSectors[i]=kFALSE;
 766   TBranch * branch=0;
 767   if (TreeH() == 0x0)
 768    {
 769      AliFatal("Can not find TreeH in folder");
 770      return;
 771    }
 772   if (fHitType>1) branch = TreeH()->GetBranch("TPC2");
 773   else branch = TreeH()->GetBranch("TPC");
 774   Stat_t ntracks = TreeH()->GetEntries();
 775   // loop over all hits
 776   AliDebug(1,Form("Got %d tracks",ntracks));
 777
 778   for(Int_t track=0;track<ntracks;track++) {
 779     ResetHits();
 780     //
 781     if (fTrackHits && fHitType&4) {
 782       TBranch * br1 = TreeH()->GetBranch("fVolumes");
 783       TBranch * br2 = TreeH()->GetBranch("fNVolumes");
 784       br1->GetEvent(track);
 785       br2->GetEvent(track);
 786       Int_t *volumes = fTrackHits->GetVolumes();
 787       for (Int_t j=0;j<fTrackHits->GetNVolumes(); j++)
 788         fActiveSectors[volumes[j]]=kTRUE;
 789     }
 790
 791     //
 792     if (fTrackHitsOld && fHitType&2) {
 793       TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
 794       br->GetEvent(track);
 795       AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
 796       for (UInt_t j=0;j<ar->GetSize();j++){
 797         fActiveSectors[((AliTrackHitsInfo*)ar->At(j))->fVolumeID] =kTRUE;
 798       }
 799     }
 800   }
 801 }
 802
 803
 804
 805
 806 //_____________________________________________________________________________
 807 void AliTPC::Digits2Raw()
 808 {
 809 // convert digits of the current event to raw data
 810
 811   static const Int_t kThreshold = 0;
 812   static const Bool_t kCompress = kTRUE;
 813
 814   fLoader->LoadDigits();
 815   TTree* digits = fLoader->TreeD();
 816   if (!digits) {
 817     AliError("No digits tree");
 818     return;
 819   }
 820
 821   AliSimDigits digarr;
 822   AliSimDigits* digrow = &digarr;
 823   digits->GetBranch("Segment")->SetAddress(&digrow);
 824
 825   const char* fileName = "AliTPCDDL.dat";
 826   AliTPCBuffer* buffer  = new AliTPCBuffer(fileName);
 827   //Verbose level
 828   // 0: Silent
 829   // 1: cout messages
 830   // 2: txt files with digits
 831   //BE CAREFUL, verbose level 2 MUST be used only for debugging and
 832   //it is highly suggested to use this mode only for debugging digits files
 833   //reasonably small, because otherwise the size of the txt files can reach
 834   //quickly several MB wasting time and disk space.
 835   buffer->SetVerbose(0);
 836
 837   Int_t nEntries = Int_t(digits->GetEntries());
 838   Int_t previousSector = -1;
 839   Int_t subSector = 0;
 840   for (Int_t i = 0; i < nEntries; i++) {
 841     digits->GetEntry(i);
 842     Int_t sector, row;
 843     fTPCParam->AdjustSectorRow(digarr.GetID(), sector, row);
 844     if(previousSector != sector) {
 845       subSector = 0;
 846       previousSector = sector;
 847     }
 848
 849     if (sector < 36) { //inner sector [0;35]
 850       if (row != 30) {
 851         //the whole row is written into the output file
 852         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
 853                                sector, subSector, row);
 854       } else {
 855         //only the pads in the range [37;48] are written into the output file
 856         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 1,
 857                                sector, subSector, row);
 858         subSector = 1;
 859         //only the pads outside the range [37;48] are written into the output file
 860         buffer->WriteRowBinary(kThreshold, digrow, 37, 48, 2,
 861                                sector, subSector, row);
 862       }//end else
 863
 864     } else { //outer sector [36;71]
 865       if (row == 54) subSector = 2;
 866       if ((row != 27) && (row != 76)) {
 867         buffer->WriteRowBinary(kThreshold, digrow, 0, 0, 0,
 868                                sector, subSector, row);
 869       } else if (row == 27) {
 870         //only the pads outside the range [43;46] are written into the output file
 871         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 2,
 872                                  sector, subSector, row);
 873         subSector = 1;
 874         //only the pads in the range [43;46] are written into the output file
 875         buffer->WriteRowBinary(kThreshold, digrow, 43, 46, 1,
 876                                  sector, subSector, row);
 877       } else if (row == 76) {
 878         //only the pads outside the range [33;88] are written into the output file
 879         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 2,
 880                                sector, subSector, row);
 881         subSector = 3;
 882         //only the pads in the range [33;88] are written into the output file
 883         buffer->WriteRowBinary(kThreshold, digrow, 33, 88, 1,
 884                                  sector, subSector, row);
 885       }
 886     }//end else
 887   }//end for
 888
 889   delete buffer;
 890   fLoader->UnloadDigits();
 891
 892   AliTPCDDLRawData rawWriter;
 893   rawWriter.SetVerbose(0);
 894
 895   rawWriter.RawData(fileName);
 896   gSystem->Unlink(fileName);
 897
 898   if (kCompress) {
 899     AliInfo("Compressing raw data");
 900     rawWriter.RawDataCompDecompress(kTRUE);
 901     gSystem->Unlink("Statistics");
 902   }
 903 }
 904
 905
 906
 907 //______________________________________________________________________
 908 AliDigitizer* AliTPC::CreateDigitizer(AliRunDigitizer* manager) const
 909 {
 910   return new AliTPCDigitizer(manager);
 911 }
 912 //__
 913 void AliTPC::SDigits2Digits2(Int_t /*eventnumber*/)
 914 {
 915   //create digits from summable digits
 916   GenerNoise(500000); //create teble with noise
 917
 918   //conect tree with sSDigits
 919   TTree *t = fLoader->TreeS();
 920
 921   if (t == 0x0) {
 922     fLoader->LoadSDigits("READ");
 923     t = fLoader->TreeS();
 924     if (t == 0x0) {
 925       AliError("Can not get input TreeS");
 926       return;
 927     }
 928   }
 929
 930   if (fLoader->TreeD() == 0x0) fLoader->MakeTree("D");
 931
 932   AliSimDigits digarr, *dummy=&digarr;
 933   TBranch* sdb = t->GetBranch("Segment");
 934   if (sdb == 0x0) {
 935     AliError("Can not find branch with segments in TreeS.");
 936     return;
 937   }
 938
 939   sdb->SetAddress(&dummy);
 940
 941   Stat_t nentries = t->GetEntries();
 942
 943   // set zero suppression
 944
 945   fTPCParam->SetZeroSup(2);
 946
 947   // get zero suppression
 948
 949   Int_t zerosup = fTPCParam->GetZeroSup();
 950
 951   //make tree with digits
 952
 953   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
 954   arr->SetClass("AliSimDigits");
 955   arr->Setup(fTPCParam);
 956   arr->MakeTree(fLoader->TreeD());
 957
 958   AliTPCParam * par = fTPCParam;
 959
 960   //Loop over segments of the TPC
 961
 962   for (Int_t n=0; n<nentries; n++) {
 963     t->GetEvent(n);
 964     Int_t sec, row;
 965     if (!par->AdjustSectorRow(digarr.GetID(),sec,row)) {
 966       AliWarning(Form("Invalid segment ID ! %d",digarr.GetID()));
 967       continue;
 968     }
 969     if (!IsSectorActive(sec)) continue;
 970
 971     AliSimDigits * digrow =(AliSimDigits*) arr->CreateRow(sec,row);
 972     Int_t nrows = digrow->GetNRows();
 973     Int_t ncols = digrow->GetNCols();
 974
 975     digrow->ExpandBuffer();
 976     digarr.ExpandBuffer();
 977     digrow->ExpandTrackBuffer();
 978     digarr.ExpandTrackBuffer();
 979
 980
 981     Short_t * pamp0 = digarr.GetDigits();
 982     Int_t   * ptracks0 = digarr.GetTracks();
 983     Short_t * pamp1 = digrow->GetDigits();
 984     Int_t   * ptracks1 = digrow->GetTracks();
 985     Int_t  nelems =nrows*ncols;
 986     Int_t saturation = fTPCParam->GetADCSat();
 987     //use internal structure of the AliDigits - for speed reason
 988     //if you cahnge implementation
 989     //of the Alidigits - it must be rewriten -
 990     for (Int_t i= 0; i<nelems; i++){
 991       Float_t q = TMath::Nint(Float_t(*pamp0)/16.+GetNoise());
 992       if (q>zerosup){
 993         if (q>saturation) q=saturation;
 994         *pamp1=(Short_t)q;
 995
 996         ptracks1[0]=ptracks0[0];
 997         ptracks1[nelems]=ptracks0[nelems];
 998         ptracks1[2*nelems]=ptracks0[2*nelems];
 999       }
1000       pamp0++;
1001       pamp1++;
1002       ptracks0++;
1003       ptracks1++;
1004     }
1005
1006     arr->StoreRow(sec,row);
1007     arr->ClearRow(sec,row);
1008   }
1009
1010
1011   //write results
1012   fLoader->WriteDigits("OVERWRITE");
1013
1014   delete arr;
1015 }
1016 //__________________________________________________________________
1017 void AliTPC::SetDefaults(){
1018   //
1019   // setting the defaults
1020   //
1021
1022   // Set response functions
1023
1024   //
1025   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1026   rl->CdGAFile();
1027   AliTPCParamSR *param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60");
1028   if(param){
1029     AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
1030     delete param;
1031     param = new AliTPCParamSR();
1032   }
1033   else {
1034     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1035   }
1036   if(!param){
1037     AliFatal("No TPC parameters found");
1038   }
1039
1040
1041   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1042   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1043   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;
1044   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1045   rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1046   rf->SetOffset(3*param->GetZSigma());
1047   rf->Update();
1048
1049   TDirectory *savedir=gDirectory;
1050   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1051   if (!f->IsOpen())
1052     AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !");
1053
1054   TString s;
1055   prfinner->Read("prf_07504_Gati_056068_d02");
1056   //PH Set different names
1057   s=prfinner->GetGRF()->GetName();
1058   s+="in";
1059   prfinner->GetGRF()->SetName(s.Data());
1060
1061   prfouter1->Read("prf_10006_Gati_047051_d03");
1062   s=prfouter1->GetGRF()->GetName();
1063   s+="out1";
1064   prfouter1->GetGRF()->SetName(s.Data());
1065
1066   prfouter2->Read("prf_15006_Gati_047051_d03");
1067   s=prfouter2->GetGRF()->GetName();
1068   s+="out2";
1069   prfouter2->GetGRF()->SetName(s.Data());
1070
1071   f->Close();
1072   savedir->cd();
1073
1074   param->SetInnerPRF(prfinner);
1075   param->SetOuter1PRF(prfouter1);
1076   param->SetOuter2PRF(prfouter2);
1077   param->SetTimeRF(rf);
1078
1079   // set fTPCParam
1080
1081   SetParam(param);
1082
1083
1084   fDefaults = 1;
1085
1086 }
1087 //__________________________________________________________________
1088 void AliTPC::Hits2Digits()
1089 {
1090   //
1091   // creates digits from hits
1092   //
1093
1094   fLoader->LoadHits("read");
1095   fLoader->LoadDigits("recreate");
1096   AliRunLoader* runLoader = fLoader->GetRunLoader();
1097
1098   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1099     runLoader->GetEvent(iEvent);
1100     SetActiveSectors();
1101     Hits2Digits(iEvent);
1102   }
1103
1104   fLoader->UnloadHits();
1105   fLoader->UnloadDigits();
1106 }
1107 //__________________________________________________________________
1108 void AliTPC::Hits2Digits(Int_t eventnumber)
1109 {
1110  //----------------------------------------------------
1111  // Loop over all sectors for a single event
1112  //----------------------------------------------------
1113   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1114   rl->GetEvent(eventnumber);
1115   if (fLoader->TreeH() == 0x0) {
1116     if(fLoader->LoadHits()) {
1117       AliError("Can not load hits.");
1118     }
1119   }
1120   SetTreeAddress();
1121
1122   if (fLoader->TreeD() == 0x0 ) {
1123     fLoader->MakeTree("D");
1124     if (fLoader->TreeD() == 0x0 ) {
1125       AliError("Can not get TreeD");
1126       return;
1127     }
1128   }
1129
1130   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1131   GenerNoise(500000); //create teble with noise
1132
1133   //setup TPCDigitsArray
1134
1135   if(GetDigitsArray()) delete GetDigitsArray();
1136
1137   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1138   arr->SetClass("AliSimDigits");
1139   arr->Setup(fTPCParam);
1140
1141   arr->MakeTree(fLoader->TreeD());
1142   SetDigitsArray(arr);
1143
1144   fDigitsSwitch=0; // standard digits
1145
1146   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1147     if (IsSectorActive(isec)) {
1148       AliDebug(1,Form("Hits2Digits","Sector %d is active.",isec));
1149       Hits2DigitsSector(isec);
1150     }
1151     else {
1152       AliDebug(1,Form("Hits2Digits","Sector %d is NOT active.",isec));
1153     }
1154
1155   fLoader->WriteDigits("OVERWRITE");
1156
1157 //this line prevents the crash in the similar one
1158 //on the beginning of this method
1159 //destructor attempts to reset the tree, which is deleted by the loader
1160 //need to be redesign
1161   if(GetDigitsArray()) delete GetDigitsArray();
1162   SetDigitsArray(0x0);
1163
1164 }
1165
1166 //__________________________________________________________________
1167 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1168 {
1169
1170   //-----------------------------------------------------------
1171   //   summable digits - 16 bit "ADC", no noise, no saturation
1172   //-----------------------------------------------------------
1173
1174   //----------------------------------------------------
1175   // Loop over all sectors for a single event
1176   //----------------------------------------------------
1177
1178   AliRunLoader* rl = fLoader->GetRunLoader();
1179
1180   rl->GetEvent(eventnumber);
1181   if (fLoader->TreeH() == 0x0) {
1182     if(fLoader->LoadHits()) {
1183       AliError("Can not load hits.");
1184       return;
1185     }
1186   }
1187   SetTreeAddress();
1188
1189
1190   if (fLoader->TreeS() == 0x0 ) {
1191     fLoader->MakeTree("S");
1192   }
1193
1194   if(fDefaults == 0) SetDefaults();
1195
1196   GenerNoise(500000); //create table with noise
1197   //setup TPCDigitsArray
1198
1199   if(GetDigitsArray()) delete GetDigitsArray();
1200
1201
1202   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1203   arr->SetClass("AliSimDigits");
1204   arr->Setup(fTPCParam);
1205   arr->MakeTree(fLoader->TreeS());
1206
1207   SetDigitsArray(arr);
1208
1209   fDigitsSwitch=1; // summable digits
1210
1211     // set zero suppression to "0"
1212
1213   fTPCParam->SetZeroSup(0);
1214
1215   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1216     if (IsSectorActive(isec)) {
1217       Hits2DigitsSector(isec);
1218     }
1219
1220   fLoader->WriteSDigits("OVERWRITE");
1221
1222 //this line prevents the crash in the similar one
1223 //on the beginning of this method
1224 //destructor attempts to reset the tree, which is deleted by the loader
1225 //need to be redesign
1226   if(GetDigitsArray()) delete GetDigitsArray();
1227   SetDigitsArray(0x0);
1228 }
1229 //__________________________________________________________________
1230
1231 void AliTPC::Hits2SDigits()
1232 {
1233
1234   //-----------------------------------------------------------
1235   //   summable digits - 16 bit "ADC", no noise, no saturation
1236   //-----------------------------------------------------------
1237
1238   fLoader->LoadHits("read");
1239   fLoader->LoadSDigits("recreate");
1240   AliRunLoader* runLoader = fLoader->GetRunLoader();
1241
1242   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1243     runLoader->GetEvent(iEvent);
1244     SetTreeAddress();
1245     SetActiveSectors();
1246     Hits2SDigits2(iEvent);
1247   }
1248
1249   fLoader->UnloadHits();
1250   fLoader->UnloadSDigits();
1251 }
1252 //_____________________________________________________________________________
1253
1254 void AliTPC::Hits2DigitsSector(Int_t isec)
1255 {
1256   //-------------------------------------------------------------------
1257   // TPC conversion from hits to digits.
1258   //-------------------------------------------------------------------
1259
1260   //-----------------------------------------------------------------
1261   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1262   //-----------------------------------------------------------------
1263
1264   //-------------------------------------------------------
1265   //  Get the access to the track hits
1266   //-------------------------------------------------------
1267
1268   // check if the parameters are set - important if one calls this method
1269   // directly, not from the Hits2Digits
1270
1271   if(fDefaults == 0) SetDefaults();
1272
1273   TTree *tH = TreeH(); // pointer to the hits tree
1274   if (tH == 0x0) {
1275     AliFatal("Can not find TreeH in folder");
1276     return;
1277   }
1278
1279   Stat_t ntracks = tH->GetEntries();
1280
1281   if( ntracks > 0){
1282
1283   //-------------------------------------------
1284   //  Only if there are any tracks...
1285   //-------------------------------------------
1286
1287     TObjArray **row;
1288
1289     Int_t nrows =fTPCParam->GetNRow(isec);
1290
1291     row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1292
1293     MakeSector(isec,nrows,tH,ntracks,row);
1294
1295     //--------------------------------------------------------
1296     //   Digitize this sector, row by row
1297     //   row[i] is the pointer to the TObjArray of TVectors,
1298     //   each one containing electrons accepted on this
1299     //   row, assigned into tracks
1300     //--------------------------------------------------------
1301
1302     Int_t i;
1303
1304     if (fDigitsArray->GetTree()==0) {
1305       AliFatal("Tree not set in fDigitsArray");
1306     }
1307
1308     for (i=0;i<nrows;i++){
1309
1310       AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1311
1312       DigitizeRow(i,isec,row);
1313
1314       fDigitsArray->StoreRow(isec,i);
1315
1316       Int_t ndig = dig->GetDigitSize();
1317
1318       AliDebug(10,
1319                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1320                     isec,i,ndig));
1321
1322       fDigitsArray->ClearRow(isec,i);
1323
1324
1325     } // end of the sector digitization
1326
1327     for(i=0;i<nrows+2;i++){
1328       row[i]->Delete();
1329       delete row[i];
1330     }
1331
1332     delete [] row; // delete the array of pointers to TObjArray-s
1333
1334   } // ntracks >0
1335
1336 } // end of Hits2DigitsSector
1337
1338
1339 //_____________________________________________________________________________
1340 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1341 {
1342   //-----------------------------------------------------------
1343   // Single row digitization, coupling from the neighbouring
1344   // rows taken into account
1345   //-----------------------------------------------------------
1346
1347   //-----------------------------------------------------------------
1348   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1349   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1350   //-----------------------------------------------------------------
1351
1352   Float_t zerosup = fTPCParam->GetZeroSup();
1353
1354   fCurrentIndex[1]= isec;
1355
1356
1357   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1358   Int_t nofTbins = fTPCParam->GetMaxTBin();
1359   Int_t indexRange[4];
1360   //
1361   //  Integrated signal for this row
1362   //  and a single track signal
1363   //
1364
1365   TMatrix *m1 = new TMatrix(0,nofPads,0,nofTbins); // integrated
1366   TMatrix *m2 = new TMatrix(0,nofPads,0,nofTbins); // single
1367   //
1368   TMatrix &total  = *m1;
1369
1370   //  Array of pointers to the label-signal list
1371
1372   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1373   Float_t  **pList = new Float_t* [nofDigits];
1374
1375   Int_t lp;
1376   Int_t i1;
1377   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1378   //
1379   //calculate signal
1380   //
1381   Int_t row1=irow;
1382   Int_t row2=irow+2;
1383   for (Int_t row= row1;row<=row2;row++){
1384     Int_t nTracks= rows[row]->GetEntries();
1385     for (i1=0;i1<nTracks;i1++){
1386       fCurrentIndex[2]= row;
1387       fCurrentIndex[3]=irow+1;
1388       if (row==irow+1){
1389         m2->Zero();  // clear single track signal matrix
1390         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1391         GetList(trackLabel,nofPads,m2,indexRange,pList);
1392       }
1393       else   GetSignal(rows[row],i1,0,m1,indexRange);
1394     }
1395   }
1396
1397   Int_t tracks[3];
1398
1399   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1400   Int_t gi=-1;
1401   Float_t fzerosup = zerosup+0.5;
1402   for(Int_t it=0;it<nofTbins;it++){
1403     for(Int_t ip=0;ip<nofPads;ip++){
1404       gi++;
1405       Float_t q=total(ip,it);
1406       if(fDigitsSwitch == 0){
1407         q+=GetNoise();
1408         if(q <=fzerosup) continue; // do not fill zeros
1409         q = TMath::Nint(q);
1410         if(q > fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat();  // saturation
1411
1412       }
1413
1414       else {
1415         if(q <= 0.) continue; // do not fill zeros
1416         if(q>2000.) q=2000.;
1417         q *= 16.;
1418         q = TMath::Nint(q);
1419       }
1420
1421       //
1422       //  "real" signal or electronic noise (list = -1)?
1423       //
1424
1425       for(Int_t j1=0;j1<3;j1++){
1426         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1427       }
1428
1429 //Begin_Html
1430 /*
1431   <A NAME="AliDigits"></A>
1432   using of AliDigits object
1433 */
1434 //End_Html
1435       dig->SetDigitFast((Short_t)q,it,ip);
1436       if (fDigitsArray->IsSimulated()) {
1437         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1438         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1439         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1440       }
1441
1442     } // end of loop over time buckets
1443   }  // end of lop over pads
1444
1445   //
1446   //  This row has been digitized, delete nonused stuff
1447   //
1448
1449   for(lp=0;lp<nofDigits;lp++){
1450     if(pList[lp]) delete [] pList[lp];
1451   }
1452
1453   delete [] pList;
1454
1455   delete m1;
1456   delete m2;
1457
1458 } // end of DigitizeRow
1459
1460 //_____________________________________________________________________________
1461
1462 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
1463              TMatrix *m1, TMatrix *m2,Int_t *indexRange)
1464 {
1465
1466   //---------------------------------------------------------------
1467   //  Calculates 2-D signal (pad,time) for a single track,
1468   //  returns a pointer to the signal matrix and the track label
1469   //  No digitization is performed at this level!!!
1470   //---------------------------------------------------------------
1471
1472   //-----------------------------------------------------------------
1473   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1474   // Modified: Marian Ivanov
1475   //-----------------------------------------------------------------
1476
1477   TVector *tv;
1478
1479   tv = (TVector*)p1->At(ntr); // pointer to a track
1480   TVector &v = *tv;
1481
1482   Float_t label = v(0);
1483   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1)-1)/2;
1484
1485   Int_t nElectrons = (tv->GetNrows()-1)/4;
1486   indexRange[0]=9999; // min pad
1487   indexRange[1]=-1; // max pad
1488   indexRange[2]=9999; //min time
1489   indexRange[3]=-1; // max time
1490
1491   TMatrix &signal = *m1;
1492   TMatrix &total = *m2;
1493   //
1494   //  Loop over all electrons
1495   //
1496   for(Int_t nel=0; nel<nElectrons; nel++){
1497     Int_t idx=nel*4;
1498     Float_t aval =  v(idx+4);
1499     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1500     Float_t xyz[3]={v(idx+1),v(idx+2),v(idx+3)};
1501     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,fCurrentIndex[3]);
1502
1503     Int_t *index = fTPCParam->GetResBin(0);
1504     Float_t *weight = & (fTPCParam->GetResWeight(0));
1505
1506     if (n>0) for (Int_t i =0; i<n; i++){
1507       Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1508
1509       if (pad>=0){
1510         Int_t time=index[2];
1511         Float_t qweight = *(weight)*eltoadcfac;
1512
1513         if (m1!=0) signal(pad,time)+=qweight;
1514         total(pad,time)+=qweight;
1515         if (indexRange[0]>pad) indexRange[0]=pad;
1516         if (indexRange[1]<pad) indexRange[1]=pad;
1517         if (indexRange[2]>time) indexRange[2]=time;
1518         if (indexRange[3]<time) indexRange[3]=time;
1519
1520         index+=3;
1521         weight++;
1522
1523       }
1524     }
1525   } // end of loop over electrons
1526
1527   return label; // returns track label when finished
1528 }
1529
1530 //_____________________________________________________________________________
1531 void AliTPC::GetList(Float_t label,Int_t np,TMatrix *m,
1532                      Int_t *indexRange, Float_t **pList)
1533 {
1534   //----------------------------------------------------------------------
1535   //  Updates the list of tracks contributing to digits for a given row
1536   //----------------------------------------------------------------------
1537
1538   //-----------------------------------------------------------------
1539   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1540   //-----------------------------------------------------------------
1541
1542   TMatrix &signal = *m;
1543
1544   // lop over nonzero digits
1545
1546   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1547     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1548
1549
1550       // accept only the contribution larger than 500 electrons (1/2 s_noise)
1551
1552       if(signal(ip,it)<0.5) continue;
1553
1554       Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1555
1556       if(!pList[globalIndex]){
1557
1558         //
1559         // Create new list (6 elements - 3 signals and 3 labels),
1560         //
1561
1562         pList[globalIndex] = new Float_t [6];
1563
1564         // set list to -1
1565
1566         *pList[globalIndex] = -1.;
1567         *(pList[globalIndex]+1) = -1.;
1568         *(pList[globalIndex]+2) = -1.;
1569         *(pList[globalIndex]+3) = -1.;
1570         *(pList[globalIndex]+4) = -1.;
1571         *(pList[globalIndex]+5) = -1.;
1572
1573         *pList[globalIndex] = label;
1574         *(pList[globalIndex]+3) = signal(ip,it);
1575       }
1576       else {
1577
1578         // check the signal magnitude
1579
1580         Float_t highest = *(pList[globalIndex]+3);
1581         Float_t middle = *(pList[globalIndex]+4);
1582         Float_t lowest = *(pList[globalIndex]+5);
1583
1584         //
1585         //  compare the new signal with already existing list
1586         //
1587
1588         if(signal(ip,it)<lowest) continue; // neglect this track
1589
1590         //
1591
1592         if (signal(ip,it)>highest){
1593           *(pList[globalIndex]+5) = middle;
1594           *(pList[globalIndex]+4) = highest;
1595           *(pList[globalIndex]+3) = signal(ip,it);
1596
1597           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1598           *(pList[globalIndex]+1) = *pList[globalIndex];
1599           *pList[globalIndex] = label;
1600         }
1601         else if (signal(ip,it)>middle){
1602           *(pList[globalIndex]+5) = middle;
1603           *(pList[globalIndex]+4) = signal(ip,it);
1604
1605           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1606           *(pList[globalIndex]+1) = label;
1607         }
1608         else{
1609           *(pList[globalIndex]+5) = signal(ip,it);
1610           *(pList[globalIndex]+2) = label;
1611         }
1612       }
1613
1614     } // end of loop over pads
1615   } // end of loop over time bins
1616
1617 }//end of GetList
1618 //___________________________________________________________________
1619 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1620                         Stat_t ntracks,TObjArray **row)
1621 {
1622
1623   //-----------------------------------------------------------------
1624   // Prepares the sector digitization, creates the vectors of
1625   // tracks for each row of this sector. The track vector
1626   // contains the track label and the position of electrons.
1627   //-----------------------------------------------------------------
1628
1629   //-----------------------------------------------------------------
1630   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1631   //-----------------------------------------------------------------
1632
1633   Float_t gasgain = fTPCParam->GetGasGain();
1634   Int_t i;
1635   Float_t xyz[4];
1636
1637   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1638   //MI change
1639   TBranch * branch=0;
1640   if (fHitType>1) branch = TH->GetBranch("TPC2");
1641   else branch = TH->GetBranch("TPC");
1642
1643
1644   //----------------------------------------------
1645   // Create TObjArray-s, one for each row,
1646   // each TObjArray will store the TVectors
1647   // of electrons, one TVectors per each track.
1648   //----------------------------------------------
1649
1650   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
1651   TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
1652
1653   for(i=0; i<nrows+2; i++){
1654     row[i] = new TObjArray;
1655     nofElectrons[i]=0;
1656     tracks[i]=0;
1657   }
1658
1659
1660
1661   //--------------------------------------------------------------------
1662   //  Loop over tracks, the "track" contains the full history
1663   //--------------------------------------------------------------------
1664
1665   Int_t previousTrack,currentTrack;
1666   previousTrack = -1; // nothing to store so far!
1667
1668   for(Int_t track=0;track<ntracks;track++){
1669     Bool_t isInSector=kTRUE;
1670     ResetHits();
1671     isInSector = TrackInVolume(isec,track);
1672     if (!isInSector) continue;
1673     //MI change
1674     branch->GetEntry(track); // get next track
1675
1676     //M.I. changes
1677
1678     tpcHit = (AliTPChit*)FirstHit(-1);
1679
1680     //--------------------------------------------------------------
1681     //  Loop over hits
1682     //--------------------------------------------------------------
1683
1684
1685     while(tpcHit){
1686
1687       Int_t sector=tpcHit->fSector; // sector number
1688       if(sector != isec){
1689         tpcHit = (AliTPChit*) NextHit();
1690         continue;
1691       }
1692
1693       currentTrack = tpcHit->Track(); // track number
1694
1695       if(currentTrack != previousTrack){
1696
1697         // store already filled fTrack
1698
1699         for(i=0;i<nrows+2;i++){
1700           if(previousTrack != -1){
1701             if(nofElectrons[i]>0){
1702               TVector &v = *tracks[i];
1703               v(0) = previousTrack;
1704               tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1705               row[i]->Add(tracks[i]);
1706             }
1707             else {
1708               delete tracks[i]; // delete empty TVector
1709               tracks[i]=0;
1710             }
1711           }
1712
1713           nofElectrons[i]=0;
1714           tracks[i] = new TVector(481); // TVectors for the next fTrack
1715
1716         } // end of loop over rows
1717
1718         previousTrack=currentTrack; // update track label
1719       }
1720
1721       Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
1722
1723       //---------------------------------------------------
1724       //  Calculate the electron attachment probability
1725       //---------------------------------------------------
1726
1727
1728       Float_t time = 1.e6*(fTPCParam->GetZLength()-TMath::Abs(tpcHit->Z()))
1729         /fTPCParam->GetDriftV();
1730       // in microseconds!
1731       Float_t attProb = fTPCParam->GetAttCoef()*
1732         fTPCParam->GetOxyCont()*time; //  fraction!
1733
1734       //-----------------------------------------------
1735       //  Loop over electrons
1736       //-----------------------------------------------
1737       Int_t index[3];
1738       index[1]=isec;
1739       for(Int_t nel=0;nel<qI;nel++){
1740         // skip if electron lost due to the attachment
1741         if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
1742         xyz[0]=tpcHit->X();
1743         xyz[1]=tpcHit->Y();
1744         xyz[2]=tpcHit->Z();
1745         //
1746         // protection for the nonphysical avalanche size (10**6 maximum)
1747         //
1748         Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
1749         xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
1750         index[0]=1;
1751
1752         TransportElectron(xyz,index);
1753         Int_t rowNumber;
1754         fTPCParam->GetPadRow(xyz,index);
1755         // row 0 - cross talk from the innermost row
1756         // row fNRow+1 cross talk from the outermost row
1757         rowNumber = index[2]+1;
1758         //transform position to local digit coordinates
1759         //relative to nearest pad row
1760         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
1761         Float_t x1,y1;
1762         if (isec <fTPCParam->GetNInnerSector()) {
1763           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
1764           y1 = fTPCParam->GetYInner(rowNumber);
1765         }
1766         else{
1767           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
1768           y1 = fTPCParam->GetYOuter(rowNumber);
1769         }
1770         // gain inefficiency at the wires edges - linear
1771         x1=TMath::Abs(x1);
1772         y1-=1.;
1773         if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.));
1774
1775         nofElectrons[rowNumber]++;
1776         //----------------------------------
1777         // Expand vector if necessary
1778         //----------------------------------
1779         if(nofElectrons[rowNumber]>120){
1780           Int_t range = tracks[rowNumber]->GetNrows();
1781           if((nofElectrons[rowNumber])>(range-1)/4){
1782
1783             tracks[rowNumber]->ResizeTo(range+400); // Add 100 electrons
1784           }
1785         }
1786
1787         TVector &v = *tracks[rowNumber];
1788         Int_t idx = 4*nofElectrons[rowNumber]-3;
1789         Real_t * position = &(((TVector&)v)(idx)); //make code faster
1790         memcpy(position,xyz,4*sizeof(Float_t));
1791
1792       } // end of loop over electrons
1793
1794       tpcHit = (AliTPChit*)NextHit();
1795
1796     } // end of loop over hits
1797   } // end of loop over tracks
1798
1799     //
1800     //   store remaining track (the last one) if not empty
1801     //
1802
1803   for(i=0;i<nrows+2;i++){
1804     if(nofElectrons[i]>0){
1805       TVector &v = *tracks[i];
1806       v(0) = previousTrack;
1807       tracks[i]->ResizeTo(4*nofElectrons[i]+1); // shrink if necessary
1808       row[i]->Add(tracks[i]);
1809     }
1810     else{
1811       delete tracks[i];
1812       tracks[i]=0;
1813     }
1814   }
1815
1816   delete [] tracks;
1817   delete [] nofElectrons;
1818
1819 } // end of MakeSector
1820
1821
1822 //_____________________________________________________________________________
1823 void AliTPC::Init()
1824 {
1825   //
1826   // Initialise TPC detector after definition of geometry
1827   //
1828   AliDebug(1,"*********************************************");
1829 }
1830
1831 //_____________________________________________________________________________
1832 void AliTPC::MakeBranch(Option_t* option)
1833 {
1834   //
1835   // Create Tree branches for the TPC.
1836   //
1837   AliDebug(1,"");
1838   Int_t buffersize = 4000;
1839   char branchname[10];
1840   sprintf(branchname,"%s",GetName());
1841
1842   const char *h = strstr(option,"H");
1843
1844   if ( h && (fHitType<=1) && (fHits == 0x0)) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
1845
1846   AliDetector::MakeBranch(option);
1847
1848   const char *d = strstr(option,"D");
1849
1850   if (fDigits   && fLoader->TreeD() && d) {
1851     MakeBranchInTree(gAlice->TreeD(), branchname, &fDigits, buffersize, 0);
1852   }
1853
1854   if (fHitType>1) MakeBranch2(option,0); // MI change 14.09.2000
1855 }
1856
1857 //_____________________________________________________________________________
1858 void AliTPC::ResetDigits()
1859 {
1860   //
1861   // Reset number of digits and the digits array for this detector
1862   //
1863   fNdigits   = 0;
1864   if (fDigits)   fDigits->Clear();
1865 }
1866
1867 //_____________________________________________________________________________
1868 void AliTPC::SetSecAL(Int_t sec)
1869 {
1870   //---------------------------------------------------
1871   // Activate/deactivate selection for lower sectors
1872   //---------------------------------------------------
1873
1874   //-----------------------------------------------------------------
1875   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1876   //-----------------------------------------------------------------
1877   fSecAL = sec;
1878 }
1879
1880 //_____________________________________________________________________________
1881 void AliTPC::SetSecAU(Int_t sec)
1882 {
1883   //----------------------------------------------------
1884   // Activate/deactivate selection for upper sectors
1885   //---------------------------------------------------
1886
1887   //-----------------------------------------------------------------
1888   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1889   //-----------------------------------------------------------------
1890   fSecAU = sec;
1891 }
1892
1893 //_____________________________________________________________________________
1894 void AliTPC::SetSecLows(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6)
1895 {
1896   //----------------------------------------
1897   // Select active lower sectors
1898   //----------------------------------------
1899
1900   //-----------------------------------------------------------------
1901   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1902   //-----------------------------------------------------------------
1903
1904   fSecLows[0] = s1;
1905   fSecLows[1] = s2;
1906   fSecLows[2] = s3;
1907   fSecLows[3] = s4;
1908   fSecLows[4] = s5;
1909   fSecLows[5] = s6;
1910 }
1911
1912 //_____________________________________________________________________________
1913 void AliTPC::SetSecUps(Int_t s1,Int_t s2,Int_t s3,Int_t s4,Int_t s5, Int_t s6,
1914                        Int_t s7, Int_t s8 ,Int_t s9 ,Int_t s10,
1915                        Int_t s11 , Int_t s12)
1916 {
1917   //--------------------------------
1918   // Select active upper sectors
1919   //--------------------------------
1920
1921   //-----------------------------------------------------------------
1922   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1923   //-----------------------------------------------------------------
1924
1925   fSecUps[0] = s1;
1926   fSecUps[1] = s2;
1927   fSecUps[2] = s3;
1928   fSecUps[3] = s4;
1929   fSecUps[4] = s5;
1930   fSecUps[5] = s6;
1931   fSecUps[6] = s7;
1932   fSecUps[7] = s8;
1933   fSecUps[8] = s9;
1934   fSecUps[9] = s10;
1935   fSecUps[10] = s11;
1936   fSecUps[11] = s12;
1937 }
1938
1939 //_____________________________________________________________________________
1940 void AliTPC::SetSens(Int_t sens)
1941 {
1942
1943   //-------------------------------------------------------------
1944   // Activates/deactivates the sensitive strips at the center of
1945   // the pad row -- this is for the space-point resolution calculations
1946   //-------------------------------------------------------------
1947
1948   //-----------------------------------------------------------------
1949   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1950   //-----------------------------------------------------------------
1951
1952   fSens = sens;
1953 }
1954
1955
1956 void AliTPC::SetSide(Float_t side=0.)
1957 {
1958   // choice of the TPC side
1959
1960   fSide = side;
1961
1962 }
1963 //____________________________________________________________________________
1964 void AliTPC::SetGasMixt(Int_t nc,Int_t c1,Int_t c2,Int_t c3,Float_t p1,
1965                            Float_t p2,Float_t p3)
1966 {
1967
1968   // gax mixture definition
1969
1970   fNoComp = nc;
1971
1972   fMixtComp[0]=c1;
1973   fMixtComp[1]=c2;
1974   fMixtComp[2]=c3;
1975
1976   fMixtProp[0]=p1;
1977   fMixtProp[1]=p2;
1978   fMixtProp[2]=p3;
1979 }
1980 //_____________________________________________________________________________
1981
1982 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
1983 {
1984   //
1985   // electron transport taking into account:
1986   // 1. diffusion,
1987   // 2.ExB at the wires
1988   // 3. nonisochronity
1989   //
1990   // xyz and index must be already transformed to system 1
1991   //
1992
1993   fTPCParam->Transform1to2(xyz,index);
1994
1995   //add diffusion
1996   Float_t driftl=xyz[2];
1997   if(driftl<0.01) driftl=0.01;
1998   driftl=TMath::Sqrt(driftl);
1999   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2000   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2001   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2002   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2003   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2004
2005   // ExB
2006
2007   if (fTPCParam->GetMWPCReadout()==kTRUE){
2008     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2009     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2010   }
2011   //add nonisochronity (not implemented yet)
2012 }
2013
2014 ClassImp(AliTPChit)
2015
2016 //_____________________________________________________________________________
2017 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits):
2018 AliHit(shunt,track)
2019 {
2020   //
2021   // Creates a TPC hit object
2022   //
2023   fSector     = vol[0];
2024   fPadRow     = vol[1];
2025   fX          = hits[0];
2026   fY          = hits[1];
2027   fZ          = hits[2];
2028   fQ          = hits[3];
2029 }
2030
2031 //________________________________________________________________________
2032 // Additional code because of the AliTPCTrackHitsV2
2033
2034 void AliTPC::MakeBranch2(Option_t *option,const char */*file*/)
2035 {
2036   //
2037   // Create a new branch in the current Root Tree
2038   // The branch of fHits is automatically split
2039   // MI change 14.09.2000
2040   AliDebug(1,"");
2041   if (fHitType<2) return;
2042   char branchname[10];
2043   sprintf(branchname,"%s2",GetName());
2044   //
2045   // Get the pointer to the header
2046   const char *cH = strstr(option,"H");
2047   //
2048   if (fTrackHits   && TreeH() && cH && fHitType&4) {
2049     AliDebug(1,"Making branch for Type 4 Hits");
2050     TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2051   }
2052
2053   if (fTrackHitsOld   && TreeH() && cH && fHitType&2) {
2054     AliDebug(1,"Making branch for Type 2 Hits");
2055     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2056                                                    TreeH(),fBufferSize,99);
2057     TreeH()->GetListOfBranches()->Add(branch);
2058   }
2059 }
2060
2061 void AliTPC::SetTreeAddress()
2062 {
2063   //Sets tree address for hits
2064   if (fHitType<=1) {
2065     if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2066     AliDetector::SetTreeAddress();
2067   }
2068   if (fHitType>1) SetTreeAddress2();
2069 }
2070
2071 void AliTPC::SetTreeAddress2()
2072 {
2073   //
2074   // Set branch address for the TrackHits Tree
2075   //
2076   AliDebug(1,"");
2077
2078   TBranch *branch;
2079   char branchname[20];
2080   sprintf(branchname,"%s2",GetName());
2081   //
2082   // Branch address for hit tree
2083   TTree *treeH = TreeH();
2084   if ((treeH)&&(fHitType&4)) {
2085     branch = treeH->GetBranch(branchname);
2086     if (branch) {
2087       branch->SetAddress(&fTrackHits);
2088       AliDebug(1,"fHitType&4 Setting");
2089     }
2090     else
2091       AliDebug(1,"fHitType&4 Failed (can not find branch)");
2092
2093   }
2094   if ((treeH)&&(fHitType&2)) {
2095     branch = treeH->GetBranch(branchname);
2096     if (branch) {
2097       branch->SetAddress(&fTrackHitsOld);
2098       AliDebug(1,"fHitType&2 Setting");
2099     }
2100     else
2101       AliDebug(1,"fHitType&2 Failed (can not find branch)");
2102   }
2103   //set address to TREETR
2104
2105   TTree *treeTR = TreeTR();
2106   if (treeTR && fTrackReferences) {
2107     branch = treeTR->GetBranch(GetName());
2108     if (branch) branch->SetAddress(&fTrackReferences);
2109   }
2110
2111 }
2112
2113 void AliTPC::FinishPrimary()
2114 {
2115   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();
2116   if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();
2117 }
2118
2119
2120 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2121 {
2122   //
2123   // add hit to the list
2124   Int_t rtrack;
2125   if (fIshunt) {
2126     int primary = gAlice->GetMCApp()->GetPrimary(track);
2127     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2128     rtrack=primary;
2129   } else {
2130     rtrack=track;
2131     gAlice->GetMCApp()->FlagTrack(track);
2132   }
2133   if (fTrackHits && fHitType&4)
2134     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2135                              hits[1],hits[2],(Int_t)hits[3]);
2136   if (fTrackHitsOld &&fHitType&2 )
2137     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2138                                 hits[1],hits[2],(Int_t)hits[3]);
2139
2140 }
2141
2142 void AliTPC::ResetHits()
2143 {
2144   if (fHitType&1) AliDetector::ResetHits();
2145   if (fHitType>1) ResetHits2();
2146 }
2147
2148 void AliTPC::ResetHits2()
2149 {
2150   //
2151   //reset hits
2152   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2153   if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2154
2155 }
2156
2157 AliHit* AliTPC::FirstHit(Int_t track)
2158 {
2159   if (fHitType>1) return FirstHit2(track);
2160   return AliDetector::FirstHit(track);
2161 }
2162 AliHit* AliTPC::NextHit()
2163 {
2164   //
2165   // gets next hit
2166   //
2167   if (fHitType>1) return NextHit2();
2168
2169   return AliDetector::NextHit();
2170 }
2171
2172 AliHit* AliTPC::FirstHit2(Int_t track)
2173 {
2174   //
2175   // Initialise the hit iterator
2176   // Return the address of the first hit for track
2177   // If track>=0 the track is read from disk
2178   // while if track<0 the first hit of the current
2179   // track is returned
2180   //
2181   if(track>=0) {
2182     gAlice->ResetHits();
2183     TreeH()->GetEvent(track);
2184   }
2185   //
2186   if (fTrackHits && fHitType&4) {
2187     fTrackHits->First();
2188     return fTrackHits->GetHit();
2189   }
2190   if (fTrackHitsOld && fHitType&2) {
2191     fTrackHitsOld->First();
2192     return fTrackHitsOld->GetHit();
2193   }
2194
2195   else return 0;
2196 }
2197
2198 AliHit* AliTPC::NextHit2()
2199 {
2200   //
2201   //Return the next hit for the current track
2202
2203
2204   if (fTrackHitsOld && fHitType&2) {
2205     fTrackHitsOld->Next();
2206     return fTrackHitsOld->GetHit();
2207   }
2208   if (fTrackHits) {
2209     fTrackHits->Next();
2210     return fTrackHits->GetHit();
2211   }
2212   else
2213     return 0;
2214 }
2215
2216 void AliTPC::LoadPoints(Int_t)
2217 {
2218   //
2219   Int_t a = 0;
2220
2221   if(fHitType==1) AliDetector::LoadPoints(a);
2222   else LoadPoints2(a);
2223 }
2224
2225
2226 void AliTPC::RemapTrackHitIDs(Int_t *map)
2227 {
2228   //
2229   // remapping
2230   //
2231   if (!fTrackHits) return;
2232
2233   if (fTrackHitsOld && fHitType&2){
2234     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2235     for (UInt_t i=0;i<arr->GetSize();i++){
2236       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2237       info->fTrackID = map[info->fTrackID];
2238     }
2239   }
2240   if (fTrackHitsOld && fHitType&4){
2241     TClonesArray * arr = fTrackHits->GetArray();;
2242     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2243       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2244       info->fTrackID = map[info->fTrackID];
2245     }
2246   }
2247 }
2248
2249 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2250 {
2251   //return bool information - is track in given volume
2252   //load only part of the track information
2253   //return true if current track is in volume
2254   //
2255   //  return kTRUE;
2256   if (fTrackHitsOld && fHitType&2) {
2257     TBranch * br = TreeH()->GetBranch("fTrackHitsInfo");
2258     br->GetEvent(track);
2259     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2260     for (UInt_t j=0;j<ar->GetSize();j++){
2261       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2262     }
2263   }
2264
2265   if (fTrackHits && fHitType&4) {
2266     TBranch * br1 = TreeH()->GetBranch("fVolumes");
2267     TBranch * br2 = TreeH()->GetBranch("fNVolumes");
2268     br2->GetEvent(track);
2269     br1->GetEvent(track);
2270     Int_t *volumes = fTrackHits->GetVolumes();
2271     Int_t nvolumes = fTrackHits->GetNVolumes();
2272     if (!volumes && nvolumes>0) {
2273       AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes));
2274       return kFALSE;
2275     }
2276     for (Int_t j=0;j<nvolumes; j++)
2277       if (volumes[j]==id) return kTRUE;
2278   }
2279
2280   if (fHitType&1) {
2281     TBranch * br = TreeH()->GetBranch("fSector");
2282     br->GetEvent(track);
2283     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2284       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2285     }
2286   }
2287   return kFALSE;
2288
2289 }
2290
2291 //_____________________________________________________________________________
2292 void AliTPC::LoadPoints2(Int_t)
2293 {
2294   //
2295   // Store x, y, z of all hits in memory
2296   //
2297   if (fTrackHits == 0 && fTrackHitsOld==0) return;
2298   //
2299   Int_t nhits =0;
2300   if (fHitType&4) nhits = fTrackHits->GetEntriesFast();
2301   if (fHitType&2) nhits = fTrackHitsOld->GetEntriesFast();
2302
2303   if (nhits == 0) return;
2304   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2305   if (fPoints == 0) fPoints = new TObjArray(tracks);
2306   AliHit *ahit;
2307   //
2308   Int_t *ntrk=new Int_t[tracks];
2309   Int_t *limi=new Int_t[tracks];
2310   Float_t **coor=new Float_t*[tracks];
2311   for(Int_t i=0;i<tracks;i++) {
2312     ntrk[i]=0;
2313     coor[i]=0;
2314     limi[i]=0;
2315   }
2316   //
2317   AliPoints *points = 0;
2318   Float_t *fp=0;
2319   Int_t trk;
2320   Int_t chunk=nhits/4+1;
2321   //
2322   // Loop over all the hits and store their position
2323   //
2324   ahit = FirstHit2(-1);
2325   while (ahit){
2326     trk=ahit->GetTrack();
2327     if(ntrk[trk]==limi[trk]) {
2328       //
2329       // Initialise a new track
2330       fp=new Float_t[3*(limi[trk]+chunk)];
2331       if(coor[trk]) {
2332         memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2333         delete [] coor[trk];
2334       }
2335       limi[trk]+=chunk;
2336       coor[trk] = fp;
2337     } else {
2338       fp = coor[trk];
2339     }
2340     fp[3*ntrk[trk]  ] = ahit->X();
2341     fp[3*ntrk[trk]+1] = ahit->Y();
2342     fp[3*ntrk[trk]+2] = ahit->Z();
2343     ntrk[trk]++;
2344     ahit = NextHit2();
2345   }
2346
2347
2348
2349   //
2350   for(trk=0; trk<tracks; ++trk) {
2351     if(ntrk[trk]) {
2352       points = new AliPoints();
2353       points->SetMarkerColor(GetMarkerColor());
2354       points->SetMarkerSize(GetMarkerSize());
2355       points->SetDetector(this);
2356       points->SetParticle(trk);
2357       points->SetPolyMarker(ntrk[trk],coor[trk],GetMarkerStyle());
2358       fPoints->AddAt(points,trk);
2359       delete [] coor[trk];
2360       coor[trk]=0;
2361     }
2362   }
2363   delete [] coor;
2364   delete [] ntrk;
2365   delete [] limi;
2366 }
2367
2368
2369 //_____________________________________________________________________________
2370 void AliTPC::LoadPoints3(Int_t)
2371 {
2372   //
2373   // Store x, y, z of all hits in memory
2374   // - only intersection point with pad row
2375   if (fTrackHits == 0) return;
2376   //
2377   Int_t nhits = fTrackHits->GetEntriesFast();
2378   if (nhits == 0) return;
2379   Int_t tracks = gAlice->GetMCApp()->GetNtrack();
2380   if (fPoints == 0) fPoints = new TObjArray(2*tracks);
2381   fPoints->Expand(2*tracks);
2382   AliHit *ahit;
2383   //
2384   Int_t *ntrk=new Int_t[tracks];
2385   Int_t *limi=new Int_t[tracks];
2386   Float_t **coor=new Float_t*[tracks];
2387   for(Int_t i=0;i<tracks;i++) {
2388     ntrk[i]=0;
2389     coor[i]=0;
2390     limi[i]=0;
2391   }
2392   //
2393   AliPoints *points = 0;
2394   Float_t *fp=0;
2395   Int_t trk;
2396   Int_t chunk=nhits/4+1;
2397   //
2398   // Loop over all the hits and store their position
2399   //
2400   ahit = FirstHit2(-1);
2401
2402   Int_t lastrow = -1;
2403   while (ahit){
2404     trk=ahit->GetTrack();
2405     Float_t  x[3]={ahit->X(),ahit->Y(),ahit->Z()};
2406     Int_t    index[3]={1,((AliTPChit*)ahit)->fSector,0};
2407     Int_t    currentrow = fTPCParam->GetPadRow(x,index) ;
2408     if (currentrow!=lastrow){
2409       lastrow = currentrow;
2410       //later calculate intersection point
2411       if(ntrk[trk]==limi[trk]) {
2412         //
2413         // Initialise a new track
2414         fp=new Float_t[3*(limi[trk]+chunk)];
2415         if(coor[trk]) {
2416           memcpy(fp,coor[trk],sizeof(Float_t)*3*limi[trk]);
2417           delete [] coor[trk];
2418         }
2419         limi[trk]+=chunk;
2420         coor[trk] = fp;
2421       } else {
2422         fp = coor[trk];
2423       }
2424       fp[3*ntrk[trk]  ] = ahit->X();
2425       fp[3*ntrk[trk]+1] = ahit->Y();
2426       fp[3*ntrk[trk]+2] = ahit->Z();
2427       ntrk[trk]++;
2428     }
2429     ahit = NextHit2();
2430   }
2431
2432   //
2433   for(trk=0; trk<tracks; ++trk) {
2434     if(ntrk[trk]) {
2435       points = new AliPoints();
2436       points->SetMarkerColor(GetMarkerColor()+1);
2437       points->SetMarkerStyle(5);
2438       points->SetMarkerSize(0.2);
2439       points->SetDetector(this);
2440       points->SetParticle(trk);
2441       points->SetPolyMarker(ntrk[trk],coor[trk],30);
2442       fPoints->AddAt(points,tracks+trk);
2443       delete [] coor[trk];
2444       coor[trk]=0;
2445     }
2446   }
2447   delete [] coor;
2448   delete [] ntrk;
2449   delete [] limi;
2450 }
2451
2452
2453
2454 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
2455 {
2456   //Makes TPC loader
2457   fLoader = new AliTPCLoader(GetName(),topfoldername);
2458   return fLoader;
2459 }
2460
2461 ////////////////////////////////////////////////////////////////////////
2462 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
2463 //
2464 // load TPC paarmeters from a given file or create new if the object
2465 // is not found there
2466 // 12/05/2003 This method should be moved to the AliTPCLoader
2467 // and one has to decide where to store the TPC parameters
2468 // M.Kowalski
2469   char paramName[50];
2470   sprintf(paramName,"75x40_100x60_150x60");
2471   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2472   if (paramTPC) {
2473     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
2474   } else {
2475     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
2476     paramTPC = new AliTPCParamSR;
2477   }
2478   return paramTPC;
2479
2480 // the older version of parameters can be accessed with this code.
2481 // In some cases, we have old parameters saved in the file but
2482 // digits were created with new parameters, it can be distinguish
2483 // by the name of TPC TreeD. The code here is just for the case
2484 // we would need to compare with old data, uncomment it if needed.
2485 //
2486 //  char paramName[50];
2487 //  sprintf(paramName,"75x40_100x60");
2488 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2489 //  if (paramTPC) {
2490 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
2491 //  } else {
2492 //    sprintf(paramName,"75x40_100x60_150x60");
2493 //    paramTPC=(AliTPCParam*)in->Get(paramName);
2494 //    if (paramTPC) {
2495 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
2496 //    } else {
2497 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2498 //          <<endl;
2499 //      paramTPC = new AliTPCParamSR;
2500 //    }
2501 //  }
2502 //  return paramTPC;
2503
2504 }
2505
2506