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