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