TPC/AliTPC.cxx

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