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