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