TPC/AliTPC.cxx

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