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   // if(param){
1214 //     AliInfo("You are using 2 pad-length geom hits with 3 pad-lenght geom digits...");
1215 //     delete param;
1216 //     param = new AliTPCParamSR();
1217 //   }
1218 //   else {
1219 //     param=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
1220 //   }
1221   param = (AliTPCParamSR*)AliTPCcalibDB::Instance()->GetParameters();
1222   if (!param->IsGeoRead()){
1223       //
1224       // read transformation matrices for gGeoManager
1225       //
1226       param->ReadGeoMatrices();
1227     }
1228   if(!param){
1229     AliFatal("No TPC parameters found");
1230   }
1231
1232
1233   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1234   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1235   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;
1236
1237
1238   //AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1239   //rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1240   //rf->SetOffset(3*param->GetZSigma());
1241   //rf->Update();
1242   //
1243   // Use gamma 4
1244   //
1245   char  strgamma4[1000];
1246   sprintf(strgamma4,"AliTPCRF1D::Gamma4((x-0.135+%f)*%f,55,160)",3*param->GetZSigma(), 1000000000*param->GetTSample()/param->GetZWidth());
1247
1248   TF1 * fgamma4 = new TF1("fgamma4",strgamma4, -1,1);
1249   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE,1000);
1250   rf->SetParam(fgamma4,param->GetZWidth(), 1,0.2);
1251   rf->SetOffset(3*param->GetZSigma());
1252   rf->Update();
1253
1254   TDirectory *savedir=gDirectory;
1255   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1256   if (!f->IsOpen())
1257     AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !");
1258
1259   TString s;
1260   prfinner->Read("prf_07504_Gati_056068_d02");
1261   //PH Set different names
1262   s=prfinner->GetGRF()->GetName();
1263   s+="in";
1264   prfinner->GetGRF()->SetName(s.Data());
1265
1266   prfouter1->Read("prf_10006_Gati_047051_d03");
1267   s=prfouter1->GetGRF()->GetName();
1268   s+="out1";
1269   prfouter1->GetGRF()->SetName(s.Data());
1270
1271   prfouter2->Read("prf_15006_Gati_047051_d03");
1272   s=prfouter2->GetGRF()->GetName();
1273   s+="out2";
1274   prfouter2->GetGRF()->SetName(s.Data());
1275
1276   f->Close();
1277   savedir->cd();
1278
1279   param->SetInnerPRF(prfinner);
1280   param->SetOuter1PRF(prfouter1);
1281   param->SetOuter2PRF(prfouter2);
1282   param->SetTimeRF(rf);
1283
1284   // set fTPCParam
1285
1286   SetParam(param);
1287
1288
1289   fDefaults = 1;
1290
1291 }
1292 //__________________________________________________________________
1293 void AliTPC::Hits2Digits()
1294 {
1295   //
1296   // creates digits from hits
1297   //
1298   if (!fTPCParam->IsGeoRead()){
1299     //
1300     // read transformation matrices for gGeoManager
1301     //
1302     fTPCParam->ReadGeoMatrices();
1303   }
1304
1305   fLoader->LoadHits("read");
1306   fLoader->LoadDigits("recreate");
1307   AliRunLoader* runLoader = fLoader->GetRunLoader();
1308
1309   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1310     //PH    runLoader->GetEvent(iEvent);
1311     Hits2Digits(iEvent);
1312   }
1313
1314   fLoader->UnloadHits();
1315   fLoader->UnloadDigits();
1316 }
1317 //__________________________________________________________________
1318 void AliTPC::Hits2Digits(Int_t eventnumber)
1319 {
1320  //----------------------------------------------------
1321  // Loop over all sectors for a single event
1322  //----------------------------------------------------
1323   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1324   rl->GetEvent(eventnumber);
1325   SetActiveSectors();
1326   if (fLoader->TreeH() == 0x0) {
1327     if(fLoader->LoadHits()) {
1328       AliError("Can not load hits.");
1329     }
1330   }
1331   SetTreeAddress();
1332
1333   if (fLoader->TreeD() == 0x0 ) {
1334     fLoader->MakeTree("D");
1335     if (fLoader->TreeD() == 0x0 ) {
1336       AliError("Can not get TreeD");
1337       return;
1338     }
1339   }
1340
1341   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1342   GenerNoise(500000); //create teble with noise
1343
1344   //setup TPCDigitsArray
1345
1346   if(GetDigitsArray()) delete GetDigitsArray();
1347
1348   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1349   arr->SetClass("AliSimDigits");
1350   arr->Setup(fTPCParam);
1351
1352   arr->MakeTree(fLoader->TreeD());
1353   SetDigitsArray(arr);
1354
1355   fDigitsSwitch=0; // standard digits
1356   // here LHC clock phase
1357   Float_t lhcph = 0.;
1358   switch (fLHCclockPhaseSw){
1359   case 0:
1360     // no phase
1361     lhcph=0.;
1362     break;
1363   case 1:
1364     // random phase
1365     lhcph = (Int_t)(gRandom->Rndm()/0.25);
1366     break;
1367   case 2:
1368     lhcph=0.;
1369     // not implemented yet
1370     break;
1371   }
1372   // adding phase to the TreeD user info
1373   fLoader->TreeD()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
1374   //
1375   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1376     if (IsSectorActive(isec)) {
1377       AliDebug(1,Form("Hits2Digits: Sector %d is active.",isec));
1378       Hits2DigitsSector(isec);
1379     }
1380     else {
1381       AliDebug(1,Form("Hits2Digits: Sector %d is NOT active.",isec));
1382     }
1383
1384   fLoader->WriteDigits("OVERWRITE");
1385
1386 //this line prevents the crash in the similar one
1387 //on the beginning of this method
1388 //destructor attempts to reset the tree, which is deleted by the loader
1389 //need to be redesign
1390   if(GetDigitsArray()) delete GetDigitsArray();
1391   SetDigitsArray(0x0);
1392
1393 }
1394
1395 //__________________________________________________________________
1396 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1397 {
1398
1399   //-----------------------------------------------------------
1400   //   summable digits - 16 bit "ADC", no noise, no saturation
1401   //-----------------------------------------------------------
1402
1403   //----------------------------------------------------
1404   // Loop over all sectors for a single event
1405   //----------------------------------------------------
1406
1407   AliRunLoader* rl = fLoader->GetRunLoader();
1408
1409   rl->GetEvent(eventnumber);
1410   if (fLoader->TreeH() == 0x0) {
1411     if(fLoader->LoadHits()) {
1412       AliError("Can not load hits.");
1413       return;
1414     }
1415   }
1416   SetTreeAddress();
1417
1418
1419   if (fLoader->TreeS() == 0x0 ) {
1420     fLoader->MakeTree("S");
1421   }
1422
1423   if(fDefaults == 0) SetDefaults();
1424
1425   GenerNoise(500000); //create table with noise
1426   //setup TPCDigitsArray
1427
1428   if(GetDigitsArray()) delete GetDigitsArray();
1429
1430
1431   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1432   arr->SetClass("AliSimDigits");
1433   arr->Setup(fTPCParam);
1434   arr->MakeTree(fLoader->TreeS());
1435
1436   SetDigitsArray(arr);
1437
1438   fDigitsSwitch=1; // summable digits
1439
1440     // set zero suppression to "0"
1441   // here LHC clock phase
1442   Float_t lhcph = 0.;
1443   switch (fLHCclockPhaseSw){
1444   case 0:
1445     // no phase
1446     lhcph=0.;
1447     break;
1448   case 1:
1449     // random phase
1450     lhcph = (Int_t)(gRandom->Rndm()/0.25);
1451     break;
1452   case 2:
1453     lhcph=0.;
1454     // not implemented yet
1455     break;
1456   }
1457   // adding phase to the TreeS user info
1458
1459   fLoader->TreeS()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
1460
1461   fTPCParam->SetZeroSup(0);
1462
1463   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1464     if (IsSectorActive(isec)) {
1465       Hits2DigitsSector(isec);
1466     }
1467
1468   fLoader->WriteSDigits("OVERWRITE");
1469
1470 //this line prevents the crash in the similar one
1471 //on the beginning of this method
1472 //destructor attempts to reset the tree, which is deleted by the loader
1473 //need to be redesign
1474   if(GetDigitsArray()) delete GetDigitsArray();
1475   SetDigitsArray(0x0);
1476 }
1477 //__________________________________________________________________
1478
1479 void AliTPC::Hits2SDigits()
1480 {
1481
1482   //-----------------------------------------------------------
1483   //   summable digits - 16 bit "ADC", no noise, no saturation
1484   //-----------------------------------------------------------
1485
1486   if (!fTPCParam->IsGeoRead()){
1487     //
1488     // read transformation matrices for gGeoManager
1489     //
1490     fTPCParam->ReadGeoMatrices();
1491   }
1492
1493   fLoader->LoadHits("read");
1494   fLoader->LoadSDigits("recreate");
1495   AliRunLoader* runLoader = fLoader->GetRunLoader();
1496
1497   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1498     runLoader->GetEvent(iEvent);
1499     SetTreeAddress();
1500     SetActiveSectors();
1501     Hits2SDigits2(iEvent);
1502   }
1503
1504   fLoader->UnloadHits();
1505   fLoader->UnloadSDigits();
1506   if (fDebugStreamer) {
1507     delete fDebugStreamer;
1508     fDebugStreamer=0;
1509   }
1510 }
1511 //_____________________________________________________________________________
1512
1513 void AliTPC::Hits2DigitsSector(Int_t isec)
1514 {
1515   //-------------------------------------------------------------------
1516   // TPC conversion from hits to digits.
1517   //-------------------------------------------------------------------
1518
1519   //-----------------------------------------------------------------
1520   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1521   //-----------------------------------------------------------------
1522
1523   //-------------------------------------------------------
1524   //  Get the access to the track hits
1525   //-------------------------------------------------------
1526
1527   // check if the parameters are set - important if one calls this method
1528   // directly, not from the Hits2Digits
1529
1530   if(fDefaults == 0) SetDefaults();
1531
1532   TTree *tH = fLoader->TreeH(); // pointer to the hits tree
1533   if (tH == 0x0) {
1534     AliFatal("Can not find TreeH in folder");
1535     return;
1536   }
1537
1538   Stat_t ntracks = tH->GetEntries();
1539
1540
1541
1542     TObjArray **row=0x0;
1543
1544     Int_t nrows =fTPCParam->GetNRow(isec);
1545
1546     row= new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1547
1548     MakeSector(isec,nrows,tH,ntracks,row);
1549
1550     //--------------------------------------------------------
1551     //   Digitize this sector, row by row
1552     //   row[i] is the pointer to the TObjArray of TVectors,
1553     //   each one containing electrons accepted on this
1554     //   row, assigned into tracks
1555     //--------------------------------------------------------
1556
1557     Int_t i;
1558
1559     if (fDigitsArray->GetTree()==0) {
1560       AliFatal("Tree not set in fDigitsArray");
1561     }
1562
1563     for (i=0;i<nrows;i++){
1564
1565       AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1566
1567       DigitizeRow(i,isec,row);
1568
1569       fDigitsArray->StoreRow(isec,i);
1570
1571       Int_t ndig = dig->GetDigitSize();
1572
1573       AliDebug(10,
1574                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1575                     isec,i,ndig));
1576
1577       fDigitsArray->ClearRow(isec,i);
1578
1579
1580     } // end of the sector digitization
1581
1582     for(i=0;i<nrows+2;i++){
1583       row[i]->Delete();
1584       delete row[i];
1585     }
1586
1587     delete [] row; // delete the array of pointers to TObjArray-s
1588
1589
1590 } // end of Hits2DigitsSector
1591
1592
1593 //_____________________________________________________________________________
1594 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1595 {
1596   //-----------------------------------------------------------
1597   // Single row digitization, coupling from the neighbouring
1598   // rows taken into account
1599   //-----------------------------------------------------------
1600
1601   //-----------------------------------------------------------------
1602   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1603   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1604   //-----------------------------------------------------------------
1605
1606   Float_t zerosup = fTPCParam->GetZeroSup();
1607   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
1608   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1609   AliTPCCalROC * gainROC = gainTPC->GetCalROC(isec);  // pad gains per given sector
1610   AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(isec);  // noise per given sector
1611
1612
1613   fCurrentIndex[1]= isec;
1614
1615
1616   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1617   Int_t nofTbins = fTPCParam->GetMaxTBin();
1618   Int_t indexRange[4];
1619   //
1620   //  Integrated signal for this row
1621   //  and a single track signal
1622   //
1623
1624   TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated
1625   TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single
1626   //
1627   TMatrixF &total  = *m1;
1628
1629   //  Array of pointers to the label-signal list
1630
1631   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1632   Float_t  **pList = new Float_t* [nofDigits];
1633
1634   Int_t lp;
1635   Int_t i1;
1636   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1637   //
1638   //calculate signal
1639   //
1640   Int_t row1=irow;
1641   Int_t row2=irow+2;
1642   for (Int_t row= row1;row<=row2;row++){
1643     Int_t nTracks= rows[row]->GetEntries();
1644     for (i1=0;i1<nTracks;i1++){
1645       fCurrentIndex[2]= row;
1646       fCurrentIndex[3]=irow+1;
1647       if (row==irow+1){
1648         m2->Zero();  // clear single track signal matrix
1649         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1650         GetList(trackLabel,nofPads,m2,indexRange,pList);
1651       }
1652       else   GetSignal(rows[row],i1,0,m1,indexRange);
1653     }
1654   }
1655
1656   Int_t tracks[3];
1657
1658   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1659   Int_t gi=-1;
1660   Float_t fzerosup = zerosup+0.5;
1661   for(Int_t it=0;it<nofTbins;it++){
1662     for(Int_t ip=0;ip<nofPads;ip++){
1663       gi++;
1664       Float_t q=total(ip,it);
1665       if(fDigitsSwitch == 0){
1666         Float_t gain = gainROC->GetValue(irow,ip);  // get gain for given - pad-row pad
1667         Float_t noisePad = noiseROC->GetValue(irow,ip);
1668         //
1669         q*=gain;
1670         q+=GetNoise()*noisePad;
1671         if(q <=fzerosup) continue; // do not fill zeros
1672         q = TMath::Nint(q);
1673         if(q >= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1;  // saturation
1674
1675       }
1676
1677       else {
1678         if(q <= 0.) continue; // do not fill zeros
1679         if(q>2000.) q=2000.;
1680         q *= 16.;
1681         q = TMath::Nint(q);
1682       }
1683
1684       //
1685       //  "real" signal or electronic noise (list = -1)?
1686       //
1687
1688       for(Int_t j1=0;j1<3;j1++){
1689         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1690       }
1691
1692 //Begin_Html
1693 /*
1694   <A NAME="AliDigits"></A>
1695   using of AliDigits object
1696 */
1697 //End_Html
1698       dig->SetDigitFast((Short_t)q,it,ip);
1699       if (fDigitsArray->IsSimulated()) {
1700         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1701         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1702         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1703       }
1704
1705     } // end of loop over time buckets
1706   }  // end of lop over pads
1707   //
1708   // test
1709   //
1710   //
1711
1712   // glitch filters if normal simulated digits
1713   //
1714   if(!fDigitsSwitch) ((AliSimDigits*)dig)->GlitchFilter();
1715   //
1716   //  This row has been digitized, delete nonused stuff
1717   //
1718
1719   for(lp=0;lp<nofDigits;lp++){
1720     if(pList[lp]) delete [] pList[lp];
1721   }
1722
1723   delete [] pList;
1724
1725   delete m1;
1726   delete m2;
1727
1728 } // end of DigitizeRow
1729
1730 //_____________________________________________________________________________
1731
1732 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
1733              TMatrixF *m1, TMatrixF *m2,Int_t *indexRange)
1734 {
1735
1736   //---------------------------------------------------------------
1737   //  Calculates 2-D signal (pad,time) for a single track,
1738   //  returns a pointer to the signal matrix and the track label
1739   //  No digitization is performed at this level!!!
1740   //---------------------------------------------------------------
1741
1742   //-----------------------------------------------------------------
1743   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1744   // Modified: Marian Ivanov
1745   //-----------------------------------------------------------------
1746
1747   TVector *tv;
1748
1749   tv = (TVector*)p1->At(ntr); // pointer to a track
1750   TVector &v = *tv;
1751
1752   Float_t label = v(0);
1753   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1))/2;
1754
1755   Int_t nElectrons = (tv->GetNrows()-1)/5;
1756   indexRange[0]=9999; // min pad
1757   indexRange[1]=-1; // max pad
1758   indexRange[2]=9999; //min time
1759   indexRange[3]=-1; // max time
1760
1761   TMatrixF &signal = *m1;
1762   TMatrixF &total = *m2;
1763   //
1764   // Get LHC clock phase
1765   //
1766   TParameter<float> *ph;
1767   if(fDigitsSwitch){// s-digits
1768     ph = (TParameter<float>*)fLoader->TreeS()->GetUserInfo()->FindObject("lhcphase0");
1769   }
1770   else{ // normal digits
1771     ph = (TParameter<float>*)fLoader->TreeD()->GetUserInfo()->FindObject("lhcphase0");
1772   }
1773   //  Loop over all electrons
1774   //
1775   for(Int_t nel=0; nel<nElectrons; nel++){
1776     Int_t idx=nel*5;
1777     Float_t aval =  v(idx+4);
1778     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1779     Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)};
1780     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,
1781                                                             fCurrentIndex[3],ph->GetVal());
1782
1783     Int_t *index = fTPCParam->GetResBin(0);
1784     Float_t *weight = & (fTPCParam->GetResWeight(0));
1785
1786     if (n>0) for (Int_t i =0; i<n; i++){
1787       Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1788
1789       if (pad>=0){
1790         Int_t time=index[2];
1791         Float_t qweight = *(weight)*eltoadcfac;
1792
1793         if (m1!=0) signal(pad,time)+=qweight;
1794         total(pad,time)+=qweight;
1795         if (indexRange[0]>pad) indexRange[0]=pad;
1796         if (indexRange[1]<pad) indexRange[1]=pad;
1797         if (indexRange[2]>time) indexRange[2]=time;
1798         if (indexRange[3]<time) indexRange[3]=time;
1799
1800         index+=3;
1801         weight++;
1802
1803       }
1804     }
1805   } // end of loop over electrons
1806
1807   return label; // returns track label when finished
1808 }
1809
1810 //_____________________________________________________________________________
1811 void AliTPC::GetList(Float_t label,Int_t np,TMatrixF *m,
1812                      Int_t *indexRange, Float_t **pList)
1813 {
1814   //----------------------------------------------------------------------
1815   //  Updates the list of tracks contributing to digits for a given row
1816   //----------------------------------------------------------------------
1817
1818   //-----------------------------------------------------------------
1819   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1820   //-----------------------------------------------------------------
1821
1822   TMatrixF &signal = *m;
1823
1824   // lop over nonzero digits
1825
1826   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1827     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1828
1829
1830       // accept only the contribution larger than 500 electrons (1/2 s_noise)
1831
1832       if(signal(ip,it)<0.5) continue;
1833
1834       Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1835
1836       if(!pList[globalIndex]){
1837
1838         //
1839         // Create new list (6 elements - 3 signals and 3 labels),
1840         //
1841
1842         pList[globalIndex] = new Float_t [6];
1843
1844         // set list to -1
1845
1846         *pList[globalIndex] = -1.;
1847         *(pList[globalIndex]+1) = -1.;
1848         *(pList[globalIndex]+2) = -1.;
1849         *(pList[globalIndex]+3) = -1.;
1850         *(pList[globalIndex]+4) = -1.;
1851         *(pList[globalIndex]+5) = -1.;
1852
1853         *pList[globalIndex] = label;
1854         *(pList[globalIndex]+3) = signal(ip,it);
1855       }
1856       else {
1857
1858         // check the signal magnitude
1859
1860         Float_t highest = *(pList[globalIndex]+3);
1861         Float_t middle = *(pList[globalIndex]+4);
1862         Float_t lowest = *(pList[globalIndex]+5);
1863
1864         //
1865         //  compare the new signal with already existing list
1866         //
1867
1868         if(signal(ip,it)<lowest) continue; // neglect this track
1869
1870         //
1871
1872         if (signal(ip,it)>highest){
1873           *(pList[globalIndex]+5) = middle;
1874           *(pList[globalIndex]+4) = highest;
1875           *(pList[globalIndex]+3) = signal(ip,it);
1876
1877           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1878           *(pList[globalIndex]+1) = *pList[globalIndex];
1879           *pList[globalIndex] = label;
1880         }
1881         else if (signal(ip,it)>middle){
1882           *(pList[globalIndex]+5) = middle;
1883           *(pList[globalIndex]+4) = signal(ip,it);
1884
1885           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1886           *(pList[globalIndex]+1) = label;
1887         }
1888         else{
1889           *(pList[globalIndex]+5) = signal(ip,it);
1890           *(pList[globalIndex]+2) = label;
1891         }
1892       }
1893
1894     } // end of loop over pads
1895   } // end of loop over time bins
1896
1897 }//end of GetList
1898 //___________________________________________________________________
1899 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1900                         Stat_t ntracks,TObjArray **row)
1901 {
1902
1903   //-----------------------------------------------------------------
1904   // Prepares the sector digitization, creates the vectors of
1905   // tracks for each row of this sector. The track vector
1906   // contains the track label and the position of electrons.
1907   //-----------------------------------------------------------------
1908
1909   //
1910   // The trasport of the electrons through TPC drift volume
1911   //    Drift (drift velocity + velocity map(not yet implemented)))
1912   //    Application of the random processes (diffusion, gas gain)
1913   //    Systematic effects (ExB effect in drift volume + ROCs)
1914   //
1915   // Algorithm:
1916   // Loop over primary electrons:
1917   //    Creation of the secondary electrons
1918   //    Loop over electrons (primary+ secondaries)
1919   //        Global coordinate frame:
1920   //          1. Skip electrons if attached
1921   //          2. ExB effect in drift volume
1922   //             a.) Simulation   calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1923   //             b.) Reconstruction -  calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1924   //          3. Generation of gas gain (Random - Exponential distribution)
1925   //          4. TransportElectron function (diffusion)
1926   //
1927   //        5. Conversion to the local coordinate frame  pad-row, pad, timebin
1928   //        6. Apply Time0 shift - AliTPCCalPad class
1929   //            a.) Plus sign in simulation
1930   //            b.) Minus sign in reconstruction
1931   // end of loop
1932   //
1933   //-----------------------------------------------------------------
1934   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1935   // Origin: Marian Ivanov,  marian.ivanov@cern.ch
1936   //-----------------------------------------------------------------
1937   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
1938   if (gAlice){ // Set correctly the magnetic field in the ExB calculation
1939     if (!calib->GetExB()){
1940       AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField());
1941       if (field) {
1942         calib->SetExBField(field);
1943       }
1944     }
1945   }
1946
1947   Float_t gasgain = fTPCParam->GetGasGain();
1948   gasgain = gasgain/fGainFactor;
1949   Int_t i;
1950   Float_t xyz[5];
1951
1952   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1953   //MI change
1954   TBranch * branch=0;
1955   if (fHitType>1) branch = TH->GetBranch("TPC2");
1956   else branch = TH->GetBranch("TPC");
1957
1958
1959   //----------------------------------------------
1960   // Create TObjArray-s, one for each row,
1961   // each TObjArray will store the TVectors
1962   // of electrons, one TVectors per each track.
1963   //----------------------------------------------
1964
1965   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
1966   TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
1967
1968   for(i=0; i<nrows+2; i++){
1969     row[i] = new TObjArray;
1970     nofElectrons[i]=0;
1971     tracks[i]=0;
1972   }
1973
1974
1975
1976   //--------------------------------------------------------------------
1977   //  Loop over tracks, the "track" contains the full history
1978   //--------------------------------------------------------------------
1979
1980   Int_t previousTrack,currentTrack;
1981   previousTrack = -1; // nothing to store so far!
1982
1983   for(Int_t track=0;track<ntracks;track++){
1984     Bool_t isInSector=kTRUE;
1985     ResetHits();
1986     isInSector = TrackInVolume(isec,track);
1987     if (!isInSector) continue;
1988     //MI change
1989     branch->GetEntry(track); // get next track
1990
1991     //M.I. changes
1992
1993     tpcHit = (AliTPChit*)FirstHit(-1);
1994
1995     //--------------------------------------------------------------
1996     //  Loop over hits
1997     //--------------------------------------------------------------
1998
1999
2000     while(tpcHit){
2001
2002       Int_t sector=tpcHit->fSector; // sector number
2003       if(sector != isec){
2004         tpcHit = (AliTPChit*) NextHit();
2005         continue;
2006       }
2007
2008       // Remove hits which arrive before the TPC opening gate signal
2009       if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
2010           /fTPCParam->GetDriftV()+tpcHit->Time())<fTPCParam->GetGateDelay()) {
2011         tpcHit = (AliTPChit*) NextHit();
2012         continue;
2013       }
2014
2015       currentTrack = tpcHit->Track(); // track number
2016
2017       if(currentTrack != previousTrack){
2018
2019         // store already filled fTrack
2020
2021         for(i=0;i<nrows+2;i++){
2022           if(previousTrack != -1){
2023             if(nofElectrons[i]>0){
2024               TVector &v = *tracks[i];
2025               v(0) = previousTrack;
2026               tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2027               row[i]->Add(tracks[i]);
2028             }
2029             else {
2030               delete tracks[i]; // delete empty TVector
2031               tracks[i]=0;
2032             }
2033           }
2034
2035           nofElectrons[i]=0;
2036           tracks[i] = new TVector(601); // TVectors for the next fTrack
2037
2038         } // end of loop over rows
2039
2040         previousTrack=currentTrack; // update track label
2041       }
2042
2043       Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2044
2045       //---------------------------------------------------
2046       //  Calculate the electron attachment probability
2047       //---------------------------------------------------
2048
2049
2050       Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
2051         /fTPCParam->GetDriftV();
2052       // in microseconds!
2053       Float_t attProb = fTPCParam->GetAttCoef()*
2054         fTPCParam->GetOxyCont()*time; //  fraction!
2055
2056       //-----------------------------------------------
2057       //  Loop over electrons
2058       //-----------------------------------------------
2059       Int_t index[3];
2060       index[1]=isec;
2061       for(Int_t nel=0;nel<qI;nel++){
2062         // skip if electron lost due to the attachment
2063         if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
2064
2065         //
2066         // ExB effect
2067         //
2068         Double_t dxyz0[3],dxyz1[3];
2069         dxyz0[0]=tpcHit->X();
2070         dxyz0[1]=tpcHit->Y();
2071         dxyz0[2]=tpcHit->Z();
2072         if (calib->GetExB()){
2073           calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
2074         }else{
2075           AliError("Not valid ExB calibration");
2076           dxyz1[0]=tpcHit->X();
2077           dxyz1[1]=tpcHit->Y();
2078           dxyz1[2]=tpcHit->Z();
2079         }
2080         xyz[0]=dxyz1[0];
2081         xyz[1]=dxyz1[1];
2082         xyz[2]=dxyz1[2];
2083         //
2084         //
2085         //
2086         // protection for the nonphysical avalanche size (10**6 maximum)
2087         //
2088         Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2089         xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
2090         index[0]=1;
2091
2092         TransportElectron(xyz,index);
2093         Int_t rowNumber;
2094         Int_t padrow = fTPCParam->GetPadRow(xyz,index);
2095         //
2096         // Add Time0 correction due unisochronity
2097         // xyz[0] - pad row coordinate
2098         // xyz[1] - pad coordinate
2099         // xyz[2] - is in now time bin coordinate system
2100         Float_t correction =0;
2101         if (calib->GetPadTime0()){
2102           if (!calib->GetPadTime0()->GetCalROC(isec)) continue;
2103           Int_t npads = fTPCParam->GetNPads(isec,padrow);
2104           //      Int_t pad  = TMath::Nint(xyz[1]+fTPCParam->GetNPads(isec,TMath::Nint(xyz[0]))*0.5);
2105           // pad numbering from -npads/2 .. npads/2-1
2106           Int_t pad  = TMath::Nint(xyz[1]+npads/2);
2107           if (pad<0) pad=0;
2108           if (pad>=npads) pad=npads-1;
2109           correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(padrow,pad);
2110           //      printf("%d\t%d\t%d\t%f\n",isec,padrow,pad,correction);
2111           if (fDebugStreamer){
2112             (*fDebugStreamer)<<"Time0"<<
2113               "isec="<<isec<<
2114               "padrow="<<padrow<<
2115               "pad="<<pad<<
2116               "x0="<<xyz[0]<<
2117               "x1="<<xyz[1]<<
2118               "x2="<<xyz[2]<<
2119               "hit.="<<tpcHit<<
2120               "cor="<<correction<<
2121               "\n";
2122           }
2123         }
2124         xyz[2]+=correction;
2125         xyz[2]+=fTPCParam->GetNTBinsL1();    // adding Level 1 time bin offset
2126         //
2127         // Electron track time (for pileup simulation)
2128         xyz[2]+=tpcHit->Time()/fTPCParam->GetTSample(); // adding time of flight
2129         xyz[4] =0;
2130
2131         //
2132         // row 0 - cross talk from the innermost row
2133         // row fNRow+1 cross talk from the outermost row
2134         rowNumber = index[2]+1;
2135         //transform position to local digit coordinates
2136         //relative to nearest pad row
2137         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
2138         /*      Float_t x1,y1;
2139         if (isec <fTPCParam->GetNInnerSector()) {
2140           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
2141           y1 = fTPCParam->GetYInner(rowNumber);
2142         }
2143         else{
2144           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
2145           y1 = fTPCParam->GetYOuter(rowNumber);
2146         }
2147         // gain inefficiency at the wires edges - linear
2148         x1=TMath::Abs(x1);
2149         y1-=1.;
2150         if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.)); */
2151
2152         nofElectrons[rowNumber]++;
2153         //----------------------------------
2154         // Expand vector if necessary
2155         //----------------------------------
2156         if(nofElectrons[rowNumber]>120){
2157           Int_t range = tracks[rowNumber]->GetNrows();
2158           if((nofElectrons[rowNumber])>(range-1)/5){
2159
2160             tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons
2161           }
2162         }
2163
2164         TVector &v = *tracks[rowNumber];
2165         Int_t idx = 5*nofElectrons[rowNumber]-4;
2166         Real_t * position = &(((TVector&)v)(idx)); //make code faster
2167         memcpy(position,xyz,5*sizeof(Float_t));
2168
2169       } // end of loop over electrons
2170
2171       tpcHit = (AliTPChit*)NextHit();
2172
2173     } // end of loop over hits
2174   } // end of loop over tracks
2175
2176     //
2177     //   store remaining track (the last one) if not empty
2178     //
2179
2180   for(i=0;i<nrows+2;i++){
2181     if(nofElectrons[i]>0){
2182       TVector &v = *tracks[i];
2183       v(0) = previousTrack;
2184       tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2185       row[i]->Add(tracks[i]);
2186     }
2187     else{
2188       delete tracks[i];
2189       tracks[i]=0;
2190     }
2191   }
2192
2193   delete [] tracks;
2194   delete [] nofElectrons;
2195
2196 } // end of MakeSector
2197
2198
2199 //_____________________________________________________________________________
2200 void AliTPC::Init()
2201 {
2202   //
2203   // Initialise TPC detector after definition of geometry
2204   //
2205   AliDebug(1,"*********************************************");
2206 }
2207
2208 //_____________________________________________________________________________
2209 void AliTPC::ResetDigits()
2210 {
2211   //
2212   // Reset number of digits and the digits array for this detector
2213   //
2214   fNdigits   = 0;
2215   if (fDigits)   fDigits->Clear();
2216 }
2217
2218
2219
2220 //_____________________________________________________________________________
2221 void AliTPC::SetSens(Int_t sens)
2222 {
2223
2224   //-------------------------------------------------------------
2225   // Activates/deactivates the sensitive strips at the center of
2226   // the pad row -- this is for the space-point resolution calculations
2227   //-------------------------------------------------------------
2228
2229   //-----------------------------------------------------------------
2230   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2231   //-----------------------------------------------------------------
2232
2233   fSens = sens;
2234 }
2235
2236
2237 void AliTPC::SetSide(Float_t side=0.)
2238 {
2239   // choice of the TPC side
2240
2241   fSide = side;
2242
2243 }
2244 //_____________________________________________________________________________
2245
2246 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2247 {
2248   //
2249   // electron transport taking into account:
2250   // 1. diffusion,
2251   // 2.ExB at the wires
2252   // 3. nonisochronity
2253   //
2254   // xyz and index must be already transformed to system 1
2255   //
2256
2257   fTPCParam->Transform1to2(xyz,index);  // mis-alignment applied in this step
2258
2259   //add diffusion
2260   Float_t driftl=xyz[2];
2261   if(driftl<0.01) driftl=0.01;
2262   driftl=TMath::Sqrt(driftl);
2263   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2264   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2265   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2266   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2267   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2268
2269   // ExB
2270
2271   if (fTPCParam->GetMWPCReadout()==kTRUE){
2272     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2273     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2274   }
2275   //add nonisochronity (not implemented yet)
2276
2277
2278 }
2279
2280 ClassImp(AliTPChit)
2281   //______________________________________________________________________
2282   AliTPChit::AliTPChit()
2283             :AliHit(),
2284              fSector(0),
2285              fPadRow(0),
2286              fQ(0),
2287              fTime(0)
2288 {
2289   //
2290   // default
2291   //
2292
2293 }
2294 //_____________________________________________________________________________
2295 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
2296           :AliHit(shunt,track),
2297              fSector(0),
2298              fPadRow(0),
2299              fQ(0),
2300              fTime(0)
2301 {
2302   //
2303   // Creates a TPC hit object
2304   //
2305   fSector     = vol[0];
2306   fPadRow     = vol[1];
2307   fX          = hits[0];
2308   fY          = hits[1];
2309   fZ          = hits[2];
2310   fQ          = hits[3];
2311   fTime       = hits[4];
2312 }
2313
2314 //________________________________________________________________________
2315 // Additional code because of the AliTPCTrackHitsV2
2316
2317 void AliTPC::MakeBranch(Option_t *option)
2318 {
2319   //
2320   // Create a new branch in the current Root Tree
2321   // The branch of fHits is automatically split
2322   // MI change 14.09.2000
2323   AliDebug(1,"");
2324   if (fHitType<2) return;
2325   char branchname[10];
2326   sprintf(branchname,"%s2",GetName());
2327   //
2328   // Get the pointer to the header
2329   const char *cH = strstr(option,"H");
2330   //
2331   if (fTrackHits   && fLoader->TreeH() && cH && fHitType&4) {
2332     AliDebug(1,"Making branch for Type 4 Hits");
2333     fLoader->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2334   }
2335
2336 //   if (fTrackHitsOld   && fLoader->TreeH() && cH && fHitType&2) {
2337 //     AliDebug(1,"Making branch for Type 2 Hits");
2338 //     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2339 //                                                    fLoader->TreeH(),fBufferSize,99);
2340 //     fLoader->TreeH()->GetListOfBranches()->Add(branch);
2341 //   }
2342 }
2343
2344 void AliTPC::SetTreeAddress()
2345 {
2346   //Sets tree address for hits
2347   if (fHitType<=1) {
2348     if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2349     AliDetector::SetTreeAddress();
2350   }
2351   if (fHitType>1) SetTreeAddress2();
2352 }
2353
2354 void AliTPC::SetTreeAddress2()
2355 {
2356   //
2357   // Set branch address for the TrackHits Tree
2358   //
2359   AliDebug(1,"");
2360
2361   TBranch *branch;
2362   char branchname[20];
2363   sprintf(branchname,"%s2",GetName());
2364   //
2365   // Branch address for hit tree
2366   TTree *treeH = fLoader->TreeH();
2367   if ((treeH)&&(fHitType&4)) {
2368     branch = treeH->GetBranch(branchname);
2369     if (branch) {
2370       branch->SetAddress(&fTrackHits);
2371       AliDebug(1,"fHitType&4 Setting");
2372     }
2373     else
2374       AliDebug(1,"fHitType&4 Failed (can not find branch)");
2375
2376   }
2377  //  if ((treeH)&&(fHitType&2)) {
2378 //     branch = treeH->GetBranch(branchname);
2379 //     if (branch) {
2380 //       branch->SetAddress(&fTrackHitsOld);
2381 //       AliDebug(1,"fHitType&2 Setting");
2382 //     }
2383 //     else
2384 //       AliDebug(1,"fHitType&2 Failed (can not find branch)");
2385 //   }
2386 }
2387
2388 void AliTPC::FinishPrimary()
2389 {
2390   if (fTrackHits &&fHitType&4)      fTrackHits->FlushHitStack();
2391   //  if (fTrackHitsOld && fHitType&2)  fTrackHitsOld->FlushHitStack();
2392 }
2393
2394
2395 void AliTPC::AddHit2(Int_t track, Int_t *vol, Float_t *hits)
2396 {
2397   //
2398   // add hit to the list
2399
2400   Int_t rtrack;
2401   if (fIshunt) {
2402     int primary = gAlice->GetMCApp()->GetPrimary(track);
2403     gAlice->GetMCApp()->Particle(primary)->SetBit(kKeepBit);
2404     rtrack=primary;
2405   } else {
2406     rtrack=track;
2407     gAlice->GetMCApp()->FlagTrack(track);
2408   }
2409   if (fTrackHits && fHitType&4)
2410     fTrackHits->AddHitKartez(vol[0],rtrack, hits[0],
2411                              hits[1],hits[2],(Int_t)hits[3],hits[4]);
2412  //  if (fTrackHitsOld &&fHitType&2 )
2413 //     fTrackHitsOld->AddHitKartez(vol[0],rtrack, hits[0],
2414 //                                 hits[1],hits[2],(Int_t)hits[3]);
2415
2416 }
2417
2418 void AliTPC::ResetHits()
2419 {
2420   if (fHitType&1) AliDetector::ResetHits();
2421   if (fHitType>1) ResetHits2();
2422 }
2423
2424 void AliTPC::ResetHits2()
2425 {
2426   //
2427   //reset hits
2428   if (fTrackHits && fHitType&4) fTrackHits->Clear();
2429   // if (fTrackHitsOld && fHitType&2) fTrackHitsOld->Clear();
2430
2431 }
2432
2433 AliHit* AliTPC::FirstHit(Int_t track)
2434 {
2435   if (fHitType>1) return FirstHit2(track);
2436   return AliDetector::FirstHit(track);
2437 }
2438 AliHit* AliTPC::NextHit()
2439 {
2440   //
2441   // gets next hit
2442   //
2443   if (fHitType>1) return NextHit2();
2444
2445   return AliDetector::NextHit();
2446 }
2447
2448 AliHit* AliTPC::FirstHit2(Int_t track)
2449 {
2450   //
2451   // Initialise the hit iterator
2452   // Return the address of the first hit for track
2453   // If track>=0 the track is read from disk
2454   // while if track<0 the first hit of the current
2455   // track is returned
2456   //
2457   if(track>=0) {
2458     gAlice->GetMCApp()->ResetHits();
2459     fLoader->TreeH()->GetEvent(track);
2460   }
2461   //
2462   if (fTrackHits && fHitType&4) {
2463     fTrackHits->First();
2464     return fTrackHits->GetHit();
2465   }
2466  //  if (fTrackHitsOld && fHitType&2) {
2467 //     fTrackHitsOld->First();
2468 //     return fTrackHitsOld->GetHit();
2469 //   }
2470
2471   else return 0;
2472 }
2473
2474 AliHit* AliTPC::NextHit2()
2475 {
2476   //
2477   //Return the next hit for the current track
2478
2479
2480 //   if (fTrackHitsOld && fHitType&2) {
2481 //     fTrackHitsOld->Next();
2482 //     return fTrackHitsOld->GetHit();
2483 //   }
2484   if (fTrackHits) {
2485     fTrackHits->Next();
2486     return fTrackHits->GetHit();
2487   }
2488   else
2489     return 0;
2490 }
2491
2492 void AliTPC::RemapTrackHitIDs(Int_t *map)
2493 {
2494   //
2495   // remapping
2496   //
2497   if (!fTrackHits) return;
2498
2499 //   if (fTrackHitsOld && fHitType&2){
2500 //     AliObjectArray * arr = fTrackHitsOld->fTrackHitsInfo;
2501 //     for (UInt_t i=0;i<arr->GetSize();i++){
2502 //       AliTrackHitsInfo * info = (AliTrackHitsInfo *)(arr->At(i));
2503 //       info->fTrackID = map[info->fTrackID];
2504 //     }
2505 //   }
2506 //  if (fTrackHitsOld && fHitType&4){
2507   if (fTrackHits && fHitType&4){
2508     TClonesArray * arr = fTrackHits->GetArray();;
2509     for (Int_t i=0;i<arr->GetEntriesFast();i++){
2510       AliTrackHitsParamV2 * info = (AliTrackHitsParamV2 *)(arr->At(i));
2511       info->SetTrackID(map[info->GetTrackID()]);
2512     }
2513   }
2514 }
2515
2516 Bool_t   AliTPC::TrackInVolume(Int_t id,Int_t track)
2517 {
2518   //return bool information - is track in given volume
2519   //load only part of the track information
2520   //return true if current track is in volume
2521   //
2522   //  return kTRUE;
2523  //  if (fTrackHitsOld && fHitType&2) {
2524 //     TBranch * br = fLoader->TreeH()->GetBranch("fTrackHitsInfo");
2525 //     br->GetEvent(track);
2526 //     AliObjectArray * ar = fTrackHitsOld->fTrackHitsInfo;
2527 //     for (UInt_t j=0;j<ar->GetSize();j++){
2528 //       if (  ((AliTrackHitsInfo*)ar->At(j))->fVolumeID==id) return kTRUE;
2529 //     }
2530 //   }
2531
2532   if (fTrackHits && fHitType&4) {
2533     TBranch * br1 = fLoader->TreeH()->GetBranch("fVolumes");
2534     TBranch * br2 = fLoader->TreeH()->GetBranch("fNVolumes");
2535     br2->GetEvent(track);
2536     br1->GetEvent(track);
2537     Int_t *volumes = fTrackHits->GetVolumes();
2538     Int_t nvolumes = fTrackHits->GetNVolumes();
2539     if (!volumes && nvolumes>0) {
2540       AliWarning(Form("Problematic track\t%d\t%d",track,nvolumes));
2541       return kFALSE;
2542     }
2543     for (Int_t j=0;j<nvolumes; j++)
2544       if (volumes[j]==id) return kTRUE;
2545   }
2546
2547   if (fHitType&1) {
2548     TBranch * br = fLoader->TreeH()->GetBranch("fSector");
2549     br->GetEvent(track);
2550     for (Int_t j=0;j<fHits->GetEntriesFast();j++){
2551       if (  ((AliTPChit*)fHits->At(j))->fSector==id) return kTRUE;
2552     }
2553   }
2554   return kFALSE;
2555
2556 }
2557
2558
2559 AliLoader* AliTPC::MakeLoader(const char* topfoldername)
2560 {
2561   //Makes TPC loader
2562   fLoader = new AliTPCLoader(GetName(),topfoldername);
2563   return fLoader;
2564 }
2565
2566 ////////////////////////////////////////////////////////////////////////
2567 AliTPCParam* AliTPC::LoadTPCParam(TFile *file) {
2568 //
2569 // load TPC paarmeters from a given file or create new if the object
2570 // is not found there
2571 // 12/05/2003 This method should be moved to the AliTPCLoader
2572 // and one has to decide where to store the TPC parameters
2573 // M.Kowalski
2574   char paramName[50];
2575   sprintf(paramName,"75x40_100x60_150x60");
2576   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2577   if (paramTPC) {
2578     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
2579   } else {
2580     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
2581     //paramTPC = new AliTPCParamSR;
2582     paramTPC = AliTPCcalibDB::Instance()->GetParameters();
2583     if (!paramTPC->IsGeoRead()){
2584       //
2585       // read transformation matrices for gGeoManager
2586       //
2587       paramTPC->ReadGeoMatrices();
2588     }
2589
2590   }
2591   return paramTPC;
2592
2593 // the older version of parameters can be accessed with this code.
2594 // In some cases, we have old parameters saved in the file but
2595 // digits were created with new parameters, it can be distinguish
2596 // by the name of TPC TreeD. The code here is just for the case
2597 // we would need to compare with old data, uncomment it if needed.
2598 //
2599 //  char paramName[50];
2600 //  sprintf(paramName,"75x40_100x60");
2601 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2602 //  if (paramTPC) {
2603 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
2604 //  } else {
2605 //    sprintf(paramName,"75x40_100x60_150x60");
2606 //    paramTPC=(AliTPCParam*)in->Get(paramName);
2607 //    if (paramTPC) {
2608 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
2609 //    } else {
2610 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2611 //          <<endl;
2612 //      paramTPC = new AliTPCParamSR;
2613 //    }
2614 //  }
2615 //  return paramTPC;
2616
2617 }
2618
2619