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