TPC/AliTPC.cxx

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