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     return;
1220   }
1221   if (!param->IsGeoRead()){
1222       //
1223       // read transformation matrices for gGeoManager
1224       //
1225       param->ReadGeoMatrices();
1226     }
1227
1228
1229
1230   AliTPCPRF2D    * prfinner   = new AliTPCPRF2D;
1231   AliTPCPRF2D    * prfouter1   = new AliTPCPRF2D;
1232   AliTPCPRF2D    * prfouter2   = new AliTPCPRF2D;
1233
1234
1235   //AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE);
1236   //rf->SetGauss(param->GetZSigma(),param->GetZWidth(),1.);
1237   //rf->SetOffset(3*param->GetZSigma());
1238   //rf->Update();
1239   //
1240   // Use gamma 4
1241   //
1242   char  strgamma4[1000];
1243   //sprintf(strgamma4,"AliTPCRF1D::Gamma4((x-0.135+%f)*%f,55,160)",3*param->GetZSigma(), 1000000000*param->GetTSample()/param->GetZWidth());
1244
1245   snprintf(strgamma4,1000,"AliTPCRF1D::Gamma4((x-0.135+%f)*%f,55,160)",3*param->GetZSigma(), 1000000000*param->GetTSample()/param->GetZWidth());
1246   TF1 * fgamma4 = new TF1("fgamma4",strgamma4, -1,1);
1247   AliTPCRF1D     * rf    = new AliTPCRF1D(kTRUE,1000);
1248   rf->SetParam(fgamma4,param->GetZWidth(), 1,0.2);
1249   rf->SetOffset(3*param->GetZSigma());
1250   rf->Update();
1251
1252   TDirectory *savedir=gDirectory;
1253   TFile *f=TFile::Open("$ALICE_ROOT/TPC/AliTPCprf2d.root");
1254   if (!f->IsOpen())
1255     AliFatal("Can't open $ALICE_ROOT/TPC/AliTPCprf2d.root !");
1256
1257   TString s;
1258   prfinner->Read("prf_07504_Gati_056068_d02");
1259   //PH Set different names
1260   s=prfinner->GetGRF()->GetName();
1261   s+="in";
1262   prfinner->GetGRF()->SetName(s.Data());
1263
1264   prfouter1->Read("prf_10006_Gati_047051_d03");
1265   s=prfouter1->GetGRF()->GetName();
1266   s+="out1";
1267   prfouter1->GetGRF()->SetName(s.Data());
1268
1269   prfouter2->Read("prf_15006_Gati_047051_d03");
1270   s=prfouter2->GetGRF()->GetName();
1271   s+="out2";
1272   prfouter2->GetGRF()->SetName(s.Data());
1273
1274   f->Close();
1275   savedir->cd();
1276
1277   param->SetInnerPRF(prfinner);
1278   param->SetOuter1PRF(prfouter1);
1279   param->SetOuter2PRF(prfouter2);
1280   param->SetTimeRF(rf);
1281
1282   // set fTPCParam
1283
1284   SetParam(param);
1285
1286
1287   fDefaults = 1;
1288
1289 }
1290 //__________________________________________________________________
1291 void AliTPC::Hits2Digits()
1292 {
1293   //
1294   // creates digits from hits
1295   //
1296   if (!fTPCParam->IsGeoRead()){
1297     //
1298     // read transformation matrices for gGeoManager
1299     //
1300     fTPCParam->ReadGeoMatrices();
1301   }
1302
1303   fLoader->LoadHits("read");
1304   fLoader->LoadDigits("recreate");
1305   AliRunLoader* runLoader = fLoader->GetRunLoader();
1306
1307   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1308     //PH    runLoader->GetEvent(iEvent);
1309     Hits2Digits(iEvent);
1310   }
1311
1312   fLoader->UnloadHits();
1313   fLoader->UnloadDigits();
1314 }
1315 //__________________________________________________________________
1316 void AliTPC::Hits2Digits(Int_t eventnumber)
1317 {
1318  //----------------------------------------------------
1319  // Loop over all sectors for a single event
1320  //----------------------------------------------------
1321   AliRunLoader* rl = (AliRunLoader*)fLoader->GetEventFolder()->FindObject(AliRunLoader::GetRunLoaderName());
1322   rl->GetEvent(eventnumber);
1323   SetActiveSectors();
1324   if (fLoader->TreeH() == 0x0) {
1325     if(fLoader->LoadHits()) {
1326       AliError("Can not load hits.");
1327     }
1328   }
1329   SetTreeAddress();
1330
1331   if (fLoader->TreeD() == 0x0 ) {
1332     fLoader->MakeTree("D");
1333     if (fLoader->TreeD() == 0x0 ) {
1334       AliError("Can not get TreeD");
1335       return;
1336     }
1337   }
1338
1339   if(fDefaults == 0) SetDefaults();  // check if the parameters are set
1340   GenerNoise(500000); //create teble with noise
1341
1342   //setup TPCDigitsArray
1343
1344   if(GetDigitsArray()) delete GetDigitsArray();
1345
1346   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1347   arr->SetClass("AliSimDigits");
1348   arr->Setup(fTPCParam);
1349
1350   arr->MakeTree(fLoader->TreeD());
1351   SetDigitsArray(arr);
1352
1353   fDigitsSwitch=0; // standard digits
1354   // here LHC clock phase
1355   Float_t lhcph = 0.;
1356   switch (fLHCclockPhaseSw){
1357   case 0:
1358     // no phase
1359     lhcph=0.;
1360     break;
1361   case 1:
1362     // random phase
1363     lhcph = (Int_t)(gRandom->Rndm()/0.25);
1364     break;
1365   case 2:
1366     lhcph=0.;
1367     // not implemented yet
1368     break;
1369   }
1370   // adding phase to the TreeD user info
1371   fLoader->TreeD()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
1372   //
1373   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1374     if (IsSectorActive(isec)) {
1375       AliDebug(1,Form("Hits2Digits: Sector %d is active.",isec));
1376       Hits2DigitsSector(isec);
1377     }
1378     else {
1379       AliDebug(1,Form("Hits2Digits: Sector %d is NOT active.",isec));
1380     }
1381
1382   fLoader->WriteDigits("OVERWRITE");
1383
1384 //this line prevents the crash in the similar one
1385 //on the beginning of this method
1386 //destructor attempts to reset the tree, which is deleted by the loader
1387 //need to be redesign
1388   if(GetDigitsArray()) delete GetDigitsArray();
1389   SetDigitsArray(0x0);
1390
1391 }
1392
1393 //__________________________________________________________________
1394 void AliTPC::Hits2SDigits2(Int_t eventnumber)
1395 {
1396
1397   //-----------------------------------------------------------
1398   //   summable digits - 16 bit "ADC", no noise, no saturation
1399   //-----------------------------------------------------------
1400
1401   //----------------------------------------------------
1402   // Loop over all sectors for a single event
1403   //----------------------------------------------------
1404
1405   AliRunLoader* rl = fLoader->GetRunLoader();
1406
1407   rl->GetEvent(eventnumber);
1408   if (fLoader->TreeH() == 0x0) {
1409     if(fLoader->LoadHits()) {
1410       AliError("Can not load hits.");
1411       return;
1412     }
1413   }
1414   SetTreeAddress();
1415
1416
1417   if (fLoader->TreeS() == 0x0 ) {
1418     fLoader->MakeTree("S");
1419   }
1420
1421   if(fDefaults == 0) SetDefaults();
1422
1423   GenerNoise(500000); //create table with noise
1424   //setup TPCDigitsArray
1425
1426   if(GetDigitsArray()) delete GetDigitsArray();
1427
1428
1429   AliTPCDigitsArray *arr = new AliTPCDigitsArray;
1430   arr->SetClass("AliSimDigits");
1431   arr->Setup(fTPCParam);
1432   arr->MakeTree(fLoader->TreeS());
1433
1434   SetDigitsArray(arr);
1435
1436   fDigitsSwitch=1; // summable digits
1437
1438     // set zero suppression to "0"
1439   // here LHC clock phase
1440   Float_t lhcph = 0.;
1441   switch (fLHCclockPhaseSw){
1442   case 0:
1443     // no phase
1444     lhcph=0.;
1445     break;
1446   case 1:
1447     // random phase
1448     lhcph = (Int_t)(gRandom->Rndm()/0.25);
1449     break;
1450   case 2:
1451     lhcph=0.;
1452     // not implemented yet
1453     break;
1454   }
1455   // adding phase to the TreeS user info
1456
1457   fLoader->TreeS()->GetUserInfo()->Add(new TParameter<float>("lhcphase0",lhcph));
1458
1459   fTPCParam->SetZeroSup(0);
1460
1461   for(Int_t isec=0;isec<fTPCParam->GetNSector();isec++)
1462     if (IsSectorActive(isec)) {
1463       Hits2DigitsSector(isec);
1464     }
1465
1466   fLoader->WriteSDigits("OVERWRITE");
1467
1468 //this line prevents the crash in the similar one
1469 //on the beginning of this method
1470 //destructor attempts to reset the tree, which is deleted by the loader
1471 //need to be redesign
1472   if(GetDigitsArray()) delete GetDigitsArray();
1473   SetDigitsArray(0x0);
1474 }
1475 //__________________________________________________________________
1476
1477 void AliTPC::Hits2SDigits()
1478 {
1479
1480   //-----------------------------------------------------------
1481   //   summable digits - 16 bit "ADC", no noise, no saturation
1482   //-----------------------------------------------------------
1483
1484   if (!fTPCParam->IsGeoRead()){
1485     //
1486     // read transformation matrices for gGeoManager
1487     //
1488     fTPCParam->ReadGeoMatrices();
1489   }
1490
1491   fLoader->LoadHits("read");
1492   fLoader->LoadSDigits("recreate");
1493   AliRunLoader* runLoader = fLoader->GetRunLoader();
1494
1495   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
1496     runLoader->GetEvent(iEvent);
1497     SetTreeAddress();
1498     SetActiveSectors();
1499     Hits2SDigits2(iEvent);
1500   }
1501
1502   fLoader->UnloadHits();
1503   fLoader->UnloadSDigits();
1504   if (fDebugStreamer) {
1505     delete fDebugStreamer;
1506     fDebugStreamer=0;
1507   }
1508 }
1509 //_____________________________________________________________________________
1510
1511 void AliTPC::Hits2DigitsSector(Int_t isec)
1512 {
1513   //-------------------------------------------------------------------
1514   // TPC conversion from hits to digits.
1515   //-------------------------------------------------------------------
1516
1517   //-----------------------------------------------------------------
1518   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1519   //-----------------------------------------------------------------
1520
1521   //-------------------------------------------------------
1522   //  Get the access to the track hits
1523   //-------------------------------------------------------
1524
1525   // check if the parameters are set - important if one calls this method
1526   // directly, not from the Hits2Digits
1527
1528   if(fDefaults == 0) SetDefaults();
1529
1530   TTree *tH = fLoader->TreeH(); // pointer to the hits tree
1531   if (tH == 0x0) {
1532     AliFatal("Can not find TreeH in folder");
1533     return;
1534   }
1535
1536   Stat_t ntracks = tH->GetEntries();
1537
1538     Int_t nrows =fTPCParam->GetNRow(isec);
1539
1540     TObjArray **row=new TObjArray* [nrows+2]; // 2 extra rows for cross talk
1541     for(Int_t j=0;j<nrows+2;j++) row[j]=0;
1542
1543     MakeSector(isec,nrows,tH,ntracks,row);
1544
1545     //--------------------------------------------------------
1546     //   Digitize this sector, row by row
1547     //   row[i] is the pointer to the TObjArray of TVectors,
1548     //   each one containing electrons accepted on this
1549     //   row, assigned into tracks
1550     //--------------------------------------------------------
1551
1552     Int_t i;
1553
1554     if (fDigitsArray->GetTree()==0) {
1555       AliFatal("Tree not set in fDigitsArray");
1556     }
1557
1558     for (i=0;i<nrows;i++){
1559
1560       AliDigits * dig = fDigitsArray->CreateRow(isec,i);
1561
1562       DigitizeRow(i,isec,row);
1563
1564       fDigitsArray->StoreRow(isec,i);
1565
1566       Int_t ndig = dig->GetDigitSize();
1567
1568       AliDebug(10,
1569                Form("*** Sector, row, compressed digits %d %d %d ***\n",
1570                     isec,i,ndig));
1571
1572       fDigitsArray->ClearRow(isec,i);
1573
1574
1575     } // end of the sector digitization
1576
1577     for(i=0;i<nrows+2;i++){
1578       row[i]->Delete();
1579       delete row[i];
1580     }
1581
1582     delete [] row; // delete the array of pointers to TObjArray-s
1583
1584
1585 } // end of Hits2DigitsSector
1586
1587
1588 //_____________________________________________________________________________
1589 void AliTPC::DigitizeRow(Int_t irow,Int_t isec,TObjArray **rows)
1590 {
1591   //-----------------------------------------------------------
1592   // Single row digitization, coupling from the neighbouring
1593   // rows taken into account
1594   //-----------------------------------------------------------
1595
1596   //-----------------------------------------------------------------
1597   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1598   // Modified: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
1599   //-----------------------------------------------------------------
1600
1601   Float_t zerosup = fTPCParam->GetZeroSup();
1602   AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetDedxGainFactor();
1603   AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise();
1604   AliTPCCalROC * gainROC = gainTPC->GetCalROC(isec);  // pad gains per given sector
1605   AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(isec);  // noise per given sector
1606
1607
1608   fCurrentIndex[1]= isec;
1609
1610
1611   Int_t nofPads = fTPCParam->GetNPads(isec,irow);
1612   Int_t nofTbins = fTPCParam->GetMaxTBin();
1613   Int_t indexRange[4];
1614   //
1615   //  Integrated signal for this row
1616   //  and a single track signal
1617   //
1618
1619   TMatrixF *m1 = new TMatrixF(0,nofPads,0,nofTbins); // integrated
1620   TMatrixF *m2 = new TMatrixF(0,nofPads,0,nofTbins); // single
1621   //
1622   TMatrixF &total  = *m1;
1623
1624   //  Array of pointers to the label-signal list
1625
1626   Int_t nofDigits = nofPads*nofTbins; // number of digits for this row
1627   Float_t  **pList = new Float_t* [nofDigits];
1628
1629   Int_t lp;
1630   Int_t i1;
1631   for(lp=0;lp<nofDigits;lp++)pList[lp]=0; // set all pointers to NULL
1632   //
1633   //calculate signal
1634   //
1635   Int_t row1=irow;
1636   Int_t row2=irow+2;
1637   for (Int_t row= row1;row<=row2;row++){
1638     Int_t nTracks= rows[row]->GetEntries();
1639     for (i1=0;i1<nTracks;i1++){
1640       fCurrentIndex[2]= row;
1641       fCurrentIndex[3]=irow+1;
1642       if (row==irow+1){
1643         m2->Zero();  // clear single track signal matrix
1644         Float_t trackLabel = GetSignal(rows[row],i1,m2,m1,indexRange);
1645         GetList(trackLabel,nofPads,m2,indexRange,pList);
1646       }
1647       else   GetSignal(rows[row],i1,0,m1,indexRange);
1648     }
1649   }
1650
1651   Int_t tracks[3];
1652
1653   AliDigits *dig = fDigitsArray->GetRow(isec,irow);
1654   Int_t gi=-1;
1655   Float_t fzerosup = zerosup+0.5;
1656   for(Int_t it=0;it<nofTbins;it++){
1657     for(Int_t ip=0;ip<nofPads;ip++){
1658       gi++;
1659       Float_t q=total(ip,it);
1660       if(fDigitsSwitch == 0){
1661         Float_t gain = gainROC->GetValue(irow,ip);  // get gain for given - pad-row pad
1662         Float_t noisePad = noiseROC->GetValue(irow,ip);
1663         //
1664         q*=gain;
1665         q+=GetNoise()*noisePad;
1666         if(q <=fzerosup) continue; // do not fill zeros
1667         q = TMath::Nint(q);
1668         if(q >= fTPCParam->GetADCSat()) q = fTPCParam->GetADCSat() - 1;  // saturation
1669
1670       }
1671
1672       else {
1673         if(q <= 0.) continue; // do not fill zeros
1674         if(q>2000.) q=2000.;
1675         q *= 16.;
1676         q = TMath::Nint(q);
1677       }
1678
1679       //
1680       //  "real" signal or electronic noise (list = -1)?
1681       //
1682
1683       for(Int_t j1=0;j1<3;j1++){
1684         tracks[j1] = (pList[gi]) ?(Int_t)(*(pList[gi]+j1)) : -2;
1685       }
1686
1687 //Begin_Html
1688 /*
1689   <A NAME="AliDigits"></A>
1690   using of AliDigits object
1691 */
1692 //End_Html
1693       dig->SetDigitFast((Short_t)q,it,ip);
1694       if (fDigitsArray->IsSimulated()) {
1695         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[0],it,ip,0);
1696         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[1],it,ip,1);
1697         ((AliSimDigits*)dig)->SetTrackIDFast(tracks[2],it,ip,2);
1698       }
1699
1700     } // end of loop over time buckets
1701   }  // end of lop over pads
1702   //
1703   // test
1704   //
1705   //
1706
1707   // glitch filters if normal simulated digits
1708   //
1709   if(!fDigitsSwitch) ((AliSimDigits*)dig)->GlitchFilter();
1710   //
1711   //  This row has been digitized, delete nonused stuff
1712   //
1713
1714   for(lp=0;lp<nofDigits;lp++){
1715     if(pList[lp]) delete [] pList[lp];
1716   }
1717
1718   delete [] pList;
1719
1720   delete m1;
1721   delete m2;
1722
1723 } // end of DigitizeRow
1724
1725 //_____________________________________________________________________________
1726
1727 Float_t AliTPC::GetSignal(TObjArray *p1, Int_t ntr,
1728              TMatrixF *m1, TMatrixF *m2,Int_t *indexRange)
1729 {
1730
1731   //---------------------------------------------------------------
1732   //  Calculates 2-D signal (pad,time) for a single track,
1733   //  returns a pointer to the signal matrix and the track label
1734   //  No digitization is performed at this level!!!
1735   //---------------------------------------------------------------
1736
1737   //-----------------------------------------------------------------
1738   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1739   // Modified: Marian Ivanov
1740   //-----------------------------------------------------------------
1741
1742   TVector *tv;
1743
1744   tv = (TVector*)p1->At(ntr); // pointer to a track
1745   TVector &v = *tv;
1746
1747   Float_t label = v(0);
1748   Int_t centralPad = (fTPCParam->GetNPads(fCurrentIndex[1],fCurrentIndex[3]-1))/2;
1749
1750   Int_t nElectrons = (tv->GetNrows()-1)/5;
1751   indexRange[0]=9999; // min pad
1752   indexRange[1]=-1; // max pad
1753   indexRange[2]=9999; //min time
1754   indexRange[3]=-1; // max time
1755
1756   TMatrixF &signal = *m1;
1757   TMatrixF &total = *m2;
1758   //
1759   // Get LHC clock phase
1760   //
1761   TParameter<float> *ph;
1762   if(fDigitsSwitch){// s-digits
1763     ph = (TParameter<float>*)fLoader->TreeS()->GetUserInfo()->FindObject("lhcphase0");
1764   }
1765   else{ // normal digits
1766     ph = (TParameter<float>*)fLoader->TreeD()->GetUserInfo()->FindObject("lhcphase0");
1767   }
1768   //  Loop over all electrons
1769   //
1770   for(Int_t nel=0; nel<nElectrons; nel++){
1771     Int_t idx=nel*5;
1772     Float_t aval =  v(idx+4);
1773     Float_t eltoadcfac=aval*fTPCParam->GetTotalNormFac();
1774     Float_t xyz[4]={v(idx+1),v(idx+2),v(idx+3),v(idx+5)};
1775     Int_t n = ((AliTPCParamSR*)fTPCParam)->CalcResponseFast(xyz,fCurrentIndex,
1776                                                             fCurrentIndex[3],ph->GetVal());
1777
1778     Int_t *index = fTPCParam->GetResBin(0);
1779     Float_t *weight = & (fTPCParam->GetResWeight(0));
1780
1781     if (n>0) for (Int_t i =0; i<n; i++){
1782       Int_t pad=index[1]+centralPad;  //in digit coordinates central pad has coordinate 0
1783
1784       if (pad>=0){
1785         Int_t time=index[2];
1786         Float_t qweight = *(weight)*eltoadcfac;
1787
1788         if (m1!=0) signal(pad,time)+=qweight;
1789         total(pad,time)+=qweight;
1790         if (indexRange[0]>pad) indexRange[0]=pad;
1791         if (indexRange[1]<pad) indexRange[1]=pad;
1792         if (indexRange[2]>time) indexRange[2]=time;
1793         if (indexRange[3]<time) indexRange[3]=time;
1794
1795         index+=3;
1796         weight++;
1797
1798       }
1799     }
1800   } // end of loop over electrons
1801
1802   return label; // returns track label when finished
1803 }
1804
1805 //_____________________________________________________________________________
1806 void AliTPC::GetList(Float_t label,Int_t np,TMatrixF *m,
1807                      Int_t *indexRange, Float_t **pList)
1808 {
1809   //----------------------------------------------------------------------
1810   //  Updates the list of tracks contributing to digits for a given row
1811   //----------------------------------------------------------------------
1812
1813   //-----------------------------------------------------------------
1814   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1815   //-----------------------------------------------------------------
1816
1817   TMatrixF &signal = *m;
1818
1819   // lop over nonzero digits
1820
1821   for(Int_t it=indexRange[2];it<indexRange[3]+1;it++){
1822     for(Int_t ip=indexRange[0];ip<indexRange[1]+1;ip++){
1823
1824
1825       // accept only the contribution larger than 500 electrons (1/2 s_noise)
1826
1827       if(signal(ip,it)<0.5) continue;
1828
1829       Int_t globalIndex = it*np+ip; // globalIndex starts from 0!
1830
1831       if(!pList[globalIndex]){
1832
1833         //
1834         // Create new list (6 elements - 3 signals and 3 labels),
1835         //
1836
1837         pList[globalIndex] = new Float_t [6];
1838
1839         // set list to -1
1840
1841         *pList[globalIndex] = -1.;
1842         *(pList[globalIndex]+1) = -1.;
1843         *(pList[globalIndex]+2) = -1.;
1844         *(pList[globalIndex]+3) = -1.;
1845         *(pList[globalIndex]+4) = -1.;
1846         *(pList[globalIndex]+5) = -1.;
1847
1848         *pList[globalIndex] = label;
1849         *(pList[globalIndex]+3) = signal(ip,it);
1850       }
1851       else {
1852
1853         // check the signal magnitude
1854
1855         Float_t highest = *(pList[globalIndex]+3);
1856         Float_t middle = *(pList[globalIndex]+4);
1857         Float_t lowest = *(pList[globalIndex]+5);
1858
1859         //
1860         //  compare the new signal with already existing list
1861         //
1862
1863         if(signal(ip,it)<lowest) continue; // neglect this track
1864
1865         //
1866
1867         if (signal(ip,it)>highest){
1868           *(pList[globalIndex]+5) = middle;
1869           *(pList[globalIndex]+4) = highest;
1870           *(pList[globalIndex]+3) = signal(ip,it);
1871
1872           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1873           *(pList[globalIndex]+1) = *pList[globalIndex];
1874           *pList[globalIndex] = label;
1875         }
1876         else if (signal(ip,it)>middle){
1877           *(pList[globalIndex]+5) = middle;
1878           *(pList[globalIndex]+4) = signal(ip,it);
1879
1880           *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
1881           *(pList[globalIndex]+1) = label;
1882         }
1883         else{
1884           *(pList[globalIndex]+5) = signal(ip,it);
1885           *(pList[globalIndex]+2) = label;
1886         }
1887       }
1888
1889     } // end of loop over pads
1890   } // end of loop over time bins
1891
1892 }//end of GetList
1893 //___________________________________________________________________
1894 void AliTPC::MakeSector(Int_t isec,Int_t nrows,TTree *TH,
1895                         Stat_t ntracks,TObjArray **row)
1896 {
1897
1898   //-----------------------------------------------------------------
1899   // Prepares the sector digitization, creates the vectors of
1900   // tracks for each row of this sector. The track vector
1901   // contains the track label and the position of electrons.
1902   //-----------------------------------------------------------------
1903
1904   //
1905   // The trasport of the electrons through TPC drift volume
1906   //    Drift (drift velocity + velocity map(not yet implemented)))
1907   //    Application of the random processes (diffusion, gas gain)
1908   //    Systematic effects (ExB effect in drift volume + ROCs)
1909   //
1910   // Algorithm:
1911   // Loop over primary electrons:
1912   //    Creation of the secondary electrons
1913   //    Loop over electrons (primary+ secondaries)
1914   //        Global coordinate frame:
1915   //          1. Skip electrons if attached
1916   //          2. ExB effect in drift volume
1917   //             a.) Simulation   calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1918   //             b.) Reconstruction -  calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
1919   //          3. Generation of gas gain (Random - Exponential distribution)
1920   //          4. TransportElectron function (diffusion)
1921   //
1922   //        5. Conversion to the local coordinate frame  pad-row, pad, timebin
1923   //        6. Apply Time0 shift - AliTPCCalPad class
1924   //            a.) Plus sign in simulation
1925   //            b.) Minus sign in reconstruction
1926   // end of loop
1927   //
1928   //-----------------------------------------------------------------
1929   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
1930   // Origin: Marian Ivanov,  marian.ivanov@cern.ch
1931   //-----------------------------------------------------------------
1932   AliTPCcalibDB* const calib=AliTPCcalibDB::Instance();
1933   if (gAlice){ // Set correctly the magnetic field in the ExB calculation
1934     if (!calib->GetExB()){
1935       AliMagF * field = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField());
1936       if (field) {
1937         calib->SetExBField(field);
1938       }
1939     }
1940   }
1941
1942   Float_t gasgain = fTPCParam->GetGasGain();
1943   gasgain = gasgain/fGainFactor;
1944   Int_t i;
1945   Float_t xyz[5];
1946
1947   AliTPChit *tpcHit; // pointer to a sigle TPC hit
1948   //MI change
1949   TBranch * branch=0;
1950   if (fHitType>1) branch = TH->GetBranch("TPC2");
1951   else branch = TH->GetBranch("TPC");
1952
1953
1954   //----------------------------------------------
1955   // Create TObjArray-s, one for each row,
1956   // each TObjArray will store the TVectors
1957   // of electrons, one TVectors per each track.
1958   //----------------------------------------------
1959
1960   Int_t *nofElectrons = new Int_t [nrows+2]; // electron counter for each row
1961   TVector **tracks = new TVector* [nrows+2]; //pointers to the track vectors
1962
1963   for(i=0; i<nrows+2; i++){
1964     row[i] = new TObjArray;
1965     nofElectrons[i]=0;
1966     tracks[i]=0;
1967   }
1968
1969
1970
1971   //--------------------------------------------------------------------
1972   //  Loop over tracks, the "track" contains the full history
1973   //--------------------------------------------------------------------
1974
1975   Int_t previousTrack,currentTrack;
1976   previousTrack = -1; // nothing to store so far!
1977
1978   for(Int_t track=0;track<ntracks;track++){
1979     Bool_t isInSector=kTRUE;
1980     ResetHits();
1981     isInSector = TrackInVolume(isec,track);
1982     if (!isInSector) continue;
1983     //MI change
1984     branch->GetEntry(track); // get next track
1985
1986     //M.I. changes
1987
1988     tpcHit = (AliTPChit*)FirstHit(-1);
1989
1990     //--------------------------------------------------------------
1991     //  Loop over hits
1992     //--------------------------------------------------------------
1993
1994
1995     while(tpcHit){
1996
1997       Int_t sector=tpcHit->fSector; // sector number
1998       if(sector != isec){
1999         tpcHit = (AliTPChit*) NextHit();
2000         continue;
2001       }
2002
2003       // Remove hits which arrive before the TPC opening gate signal
2004       if(((fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
2005           /fTPCParam->GetDriftV()+tpcHit->Time())<fTPCParam->GetGateDelay()) {
2006         tpcHit = (AliTPChit*) NextHit();
2007         continue;
2008       }
2009
2010       currentTrack = tpcHit->Track(); // track number
2011
2012       if(currentTrack != previousTrack){
2013
2014         // store already filled fTrack
2015
2016         for(i=0;i<nrows+2;i++){
2017           if(previousTrack != -1){
2018             if(nofElectrons[i]>0){
2019               TVector &v = *tracks[i];
2020               v(0) = previousTrack;
2021               tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2022               row[i]->Add(tracks[i]);
2023             }
2024             else {
2025               delete tracks[i]; // delete empty TVector
2026               tracks[i]=0;
2027             }
2028           }
2029
2030           nofElectrons[i]=0;
2031           tracks[i] = new TVector(601); // TVectors for the next fTrack
2032
2033         } // end of loop over rows
2034
2035         previousTrack=currentTrack; // update track label
2036       }
2037
2038       Int_t qI = (Int_t) (tpcHit->fQ); // energy loss (number of electrons)
2039
2040       //---------------------------------------------------
2041       //  Calculate the electron attachment probability
2042       //---------------------------------------------------
2043
2044
2045       Float_t time = 1.e6*(fTPCParam->GetZLength(isec)-TMath::Abs(tpcHit->Z()))
2046         /fTPCParam->GetDriftV();
2047       // in microseconds!
2048       Float_t attProb = fTPCParam->GetAttCoef()*
2049         fTPCParam->GetOxyCont()*time; //  fraction!
2050
2051       //-----------------------------------------------
2052       //  Loop over electrons
2053       //-----------------------------------------------
2054       Int_t index[3];
2055       index[1]=isec;
2056       for(Int_t nel=0;nel<qI;nel++){
2057         // skip if electron lost due to the attachment
2058         if((gRandom->Rndm(0)) < attProb) continue; // electron lost!
2059
2060         //
2061         // ExB effect
2062         //
2063         Double_t dxyz0[3],dxyz1[3];
2064         dxyz0[0]=tpcHit->X();
2065         dxyz0[1]=tpcHit->Y();
2066         dxyz0[2]=tpcHit->Z();
2067         if (calib->GetExB()){
2068           calib->GetExB()->CorrectInverse(dxyz0,dxyz1);
2069         }else{
2070           AliError("Not valid ExB calibration");
2071           dxyz1[0]=tpcHit->X();
2072           dxyz1[1]=tpcHit->Y();
2073           dxyz1[2]=tpcHit->Z();
2074         }
2075         xyz[0]=dxyz1[0];
2076         xyz[1]=dxyz1[1];
2077         xyz[2]=dxyz1[2];
2078         //
2079         //
2080         //
2081         // protection for the nonphysical avalanche size (10**6 maximum)
2082         //
2083         Double_t rn=TMath::Max(gRandom->Rndm(0),1.93e-22);
2084         xyz[3]= (Float_t) (-gasgain*TMath::Log(rn));
2085         index[0]=1;
2086
2087         TransportElectron(xyz,index);
2088         Int_t rowNumber;
2089         Int_t padrow = fTPCParam->GetPadRow(xyz,index);
2090         //
2091         // Add Time0 correction due unisochronity
2092         // xyz[0] - pad row coordinate
2093         // xyz[1] - pad coordinate
2094         // xyz[2] - is in now time bin coordinate system
2095         Float_t correction =0;
2096         if (calib->GetPadTime0()){
2097           if (!calib->GetPadTime0()->GetCalROC(isec)) continue;
2098           Int_t npads = fTPCParam->GetNPads(isec,padrow);
2099           //      Int_t pad  = TMath::Nint(xyz[1]+fTPCParam->GetNPads(isec,TMath::Nint(xyz[0]))*0.5);
2100           // pad numbering from -npads/2 .. npads/2-1
2101           Int_t pad  = TMath::Nint(xyz[1]+npads/2);
2102           if (pad<0) pad=0;
2103           if (pad>=npads) pad=npads-1;
2104           correction = calib->GetPadTime0()->GetCalROC(isec)->GetValue(padrow,pad);
2105           //      printf("%d\t%d\t%d\t%f\n",isec,padrow,pad,correction);
2106           if (fDebugStreamer){
2107             (*fDebugStreamer)<<"Time0"<<
2108               "isec="<<isec<<
2109               "padrow="<<padrow<<
2110               "pad="<<pad<<
2111               "x0="<<xyz[0]<<
2112               "x1="<<xyz[1]<<
2113               "x2="<<xyz[2]<<
2114               "hit.="<<tpcHit<<
2115               "cor="<<correction<<
2116               "\n";
2117           }
2118         }
2119         xyz[2]+=correction;
2120         xyz[2]+=fTPCParam->GetNTBinsL1();    // adding Level 1 time bin offset
2121         //
2122         // Electron track time (for pileup simulation)
2123         xyz[2]+=tpcHit->Time()/fTPCParam->GetTSample(); // adding time of flight
2124         xyz[4] =0;
2125
2126         //
2127         // row 0 - cross talk from the innermost row
2128         // row fNRow+1 cross talk from the outermost row
2129         rowNumber = index[2]+1;
2130         //transform position to local digit coordinates
2131         //relative to nearest pad row
2132         if ((rowNumber<0)||rowNumber>fTPCParam->GetNRow(isec)+1) continue;
2133         /*      Float_t x1,y1;
2134         if (isec <fTPCParam->GetNInnerSector()) {
2135           x1 = xyz[1]*fTPCParam->GetInnerPadPitchWidth();
2136           y1 = fTPCParam->GetYInner(rowNumber);
2137         }
2138         else{
2139           x1=xyz[1]*fTPCParam->GetOuterPadPitchWidth();
2140           y1 = fTPCParam->GetYOuter(rowNumber);
2141         }
2142         // gain inefficiency at the wires edges - linear
2143         x1=TMath::Abs(x1);
2144         y1-=1.;
2145         if(x1>y1) xyz[3]*=TMath::Max(1.e-6,(y1-x1+1.)); */
2146
2147         nofElectrons[rowNumber]++;
2148         //----------------------------------
2149         // Expand vector if necessary
2150         //----------------------------------
2151         if(nofElectrons[rowNumber]>120){
2152           Int_t range = tracks[rowNumber]->GetNrows();
2153           if((nofElectrons[rowNumber])>(range-1)/5){
2154
2155             tracks[rowNumber]->ResizeTo(range+500); // Add 100 electrons
2156           }
2157         }
2158
2159         TVector &v = *tracks[rowNumber];
2160         Int_t idx = 5*nofElectrons[rowNumber]-4;
2161         Real_t * position = &(((TVector&)v)(idx)); //make code faster
2162         memcpy(position,xyz,5*sizeof(Float_t));
2163
2164       } // end of loop over electrons
2165
2166       tpcHit = (AliTPChit*)NextHit();
2167
2168     } // end of loop over hits
2169   } // end of loop over tracks
2170
2171     //
2172     //   store remaining track (the last one) if not empty
2173     //
2174
2175   for(i=0;i<nrows+2;i++){
2176     if(nofElectrons[i]>0){
2177       TVector &v = *tracks[i];
2178       v(0) = previousTrack;
2179       tracks[i]->ResizeTo(5*nofElectrons[i]+1); // shrink if necessary
2180       row[i]->Add(tracks[i]);
2181     }
2182     else{
2183       delete tracks[i];
2184       tracks[i]=0;
2185     }
2186   }
2187
2188   delete [] tracks;
2189   delete [] nofElectrons;
2190
2191 } // end of MakeSector
2192
2193
2194 //_____________________________________________________________________________
2195 void AliTPC::Init()
2196 {
2197   //
2198   // Initialise TPC detector after definition of geometry
2199   //
2200   AliDebug(1,"*********************************************");
2201 }
2202
2203 //_____________________________________________________________________________
2204 void AliTPC::ResetDigits()
2205 {
2206   //
2207   // Reset number of digits and the digits array for this detector
2208   //
2209   fNdigits   = 0;
2210   if (fDigits)   fDigits->Clear();
2211 }
2212
2213
2214
2215 //_____________________________________________________________________________
2216 void AliTPC::SetSens(Int_t sens)
2217 {
2218
2219   //-------------------------------------------------------------
2220   // Activates/deactivates the sensitive strips at the center of
2221   // the pad row -- this is for the space-point resolution calculations
2222   //-------------------------------------------------------------
2223
2224   //-----------------------------------------------------------------
2225   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
2226   //-----------------------------------------------------------------
2227
2228   fSens = sens;
2229 }
2230
2231
2232 void AliTPC::SetSide(Float_t side=0.)
2233 {
2234   // choice of the TPC side
2235
2236   fSide = side;
2237
2238 }
2239 //_____________________________________________________________________________
2240
2241 void AliTPC::TransportElectron(Float_t *xyz, Int_t *index)
2242 {
2243   //
2244   // electron transport taking into account:
2245   // 1. diffusion,
2246   // 2.ExB at the wires
2247   // 3. nonisochronity
2248   //
2249   // xyz and index must be already transformed to system 1
2250   //
2251
2252   fTPCParam->Transform1to2(xyz,index);  // mis-alignment applied in this step
2253
2254   //add diffusion
2255   Float_t driftl=xyz[2];
2256   if(driftl<0.01) driftl=0.01;
2257   driftl=TMath::Sqrt(driftl);
2258   Float_t sigT = driftl*(fTPCParam->GetDiffT());
2259   Float_t sigL = driftl*(fTPCParam->GetDiffL());
2260   xyz[0]=gRandom->Gaus(xyz[0],sigT);
2261   xyz[1]=gRandom->Gaus(xyz[1],sigT);
2262   xyz[2]=gRandom->Gaus(xyz[2],sigL);
2263
2264   // ExB
2265
2266   if (fTPCParam->GetMWPCReadout()==kTRUE){
2267     Float_t dx = fTPCParam->Transform2to2NearestWire(xyz,index);
2268     xyz[1]+=dx*(fTPCParam->GetOmegaTau());
2269   }
2270   //add nonisochronity (not implemented yet)
2271
2272
2273 }
2274
2275 ClassImp(AliTPChit)
2276   //______________________________________________________________________
2277   AliTPChit::AliTPChit()
2278             :AliHit(),
2279              fSector(0),
2280              fPadRow(0),
2281              fQ(0),
2282              fTime(0)
2283 {
2284   //
2285   // default
2286   //
2287
2288 }
2289 //_____________________________________________________________________________
2290 AliTPChit::AliTPChit(Int_t shunt, Int_t track, Int_t *vol, Float_t *hits)
2291           :AliHit(shunt,track),
2292              fSector(0),
2293              fPadRow(0),
2294              fQ(0),
2295              fTime(0)
2296 {
2297   //
2298   // Creates a TPC hit object
2299   //
2300   fSector     = vol[0];
2301   fPadRow     = vol[1];
2302   fX          = hits[0];
2303   fY          = hits[1];
2304   fZ          = hits[2];
2305   fQ          = hits[3];
2306   fTime       = hits[4];
2307 }
2308
2309 //________________________________________________________________________
2310 // Additional code because of the AliTPCTrackHitsV2
2311
2312 void AliTPC::MakeBranch(Option_t *option)
2313 {
2314   //
2315   // Create a new branch in the current Root Tree
2316   // The branch of fHits is automatically split
2317   // MI change 14.09.2000
2318   AliDebug(1,"");
2319   if (fHitType<2) return;
2320   char branchname[10];
2321   //sprintf(branchname,"%s2",GetName());
2322   snprintf(branchname,10,"%s2",GetName());
2323   //
2324   // Get the pointer to the header
2325   const char *cH = strstr(option,"H");
2326   //
2327   if (fTrackHits   && fLoader->TreeH() && cH && fHitType&4) {
2328     AliDebug(1,"Making branch for Type 4 Hits");
2329     fLoader->TreeH()->Branch(branchname,"AliTPCTrackHitsV2",&fTrackHits,fBufferSize,99);
2330   }
2331
2332 //   if (fTrackHitsOld   && fLoader->TreeH() && cH && fHitType&2) {
2333 //     AliDebug(1,"Making branch for Type 2 Hits");
2334 //     AliObjectBranch * branch = new AliObjectBranch(branchname,"AliTPCTrackHits",&fTrackHitsOld,
2335 //                                                    fLoader->TreeH(),fBufferSize,99);
2336 //     fLoader->TreeH()->GetListOfBranches()->Add(branch);
2337 //   }
2338 }
2339
2340 void AliTPC::SetTreeAddress()
2341 {
2342   //Sets tree address for hits
2343   if (fHitType<=1) {
2344     if (fHits == 0x0 ) fHits = new TClonesArray("AliTPChit", 176);//skowron 20.06.03
2345     AliDetector::SetTreeAddress();
2346   }
2347   if (fHitType>1) SetTreeAddress2();
2348 }
2349
2350 void AliTPC::SetTreeAddress2()
2351 {
2352   //
2353   // Set branch address for the TrackHits Tree
2354   //
2355   AliDebug(1,"");
2356
2357   TBranch *branch;
2358   char branchname[20];
2359   //sprintf(branchname,"%s2",GetName());
2360   snprintf(branchname,20,"%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   snprintf(paramName,50,"75x40_100x60_150x60");
2574   AliTPCParam *paramTPC=(AliTPCParam*)file->Get(paramName);
2575   if (paramTPC) {
2576     AliDebugClass(1,Form("TPC parameters %s found.",paramName));
2577   } else {
2578     AliWarningClass("TPC parameters not found. Create new (they may be incorrect)");
2579     //paramTPC = new AliTPCParamSR;
2580     paramTPC = AliTPCcalibDB::Instance()->GetParameters();
2581     if (!paramTPC->IsGeoRead()){
2582       //
2583       // read transformation matrices for gGeoManager
2584       //
2585       paramTPC->ReadGeoMatrices();
2586     }
2587
2588   }
2589   return paramTPC;
2590
2591 // the older version of parameters can be accessed with this code.
2592 // In some cases, we have old parameters saved in the file but
2593 // digits were created with new parameters, it can be distinguish
2594 // by the name of TPC TreeD. The code here is just for the case
2595 // we would need to compare with old data, uncomment it if needed.
2596 //
2597 //  char paramName[50];
2598 //  sprintf(paramName,"75x40_100x60");
2599 //  AliTPCParam *paramTPC=(AliTPCParam*)in->Get(paramName);
2600 //  if (paramTPC) {
2601 //    cout<<"TPC parameters "<<paramName<<" found."<<endl;
2602 //  } else {
2603 //    sprintf(paramName,"75x40_100x60_150x60");
2604 //    paramTPC=(AliTPCParam*)in->Get(paramName);
2605 //    if (paramTPC) {
2606 //      cout<<"TPC parameters "<<paramName<<" found."<<endl;
2607 //    } else {
2608 //      cerr<<"TPC parameters not found. Create new (they may be incorrect)."
2609 //          <<endl;
2610 //      paramTPC = new AliTPCParamSR;
2611 //    }
2612 //  }
2613 //  return paramTPC;
2614
2615 }
2616
2617