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