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