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