TPC/AliTPC.cxx

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