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