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