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