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