TPC/AliTPCFast.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 //  fast TPC cluster simulation                                              //
  21 //                                                                          //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24 #include <TParticle.h>
  25 #include <TVector.h>
  26 #include <TRandom.h>
  27 #include <TTree.h>
  28
  29 #include "AliRunLoader.h"
  30 #include "AliLoader.h"
  31 #include "AliRun.h"
  32 #include "AliMC.h"
  33 #include "AliTPC.h"
  34 #include "AliTPCParam.h"
  35 #include "AliTPCClustersArray.h"
  36 #include "AliTPCClustersRow.h"
  37 #include "AliTPCcluster.h"
  38 #include "AliComplexCluster.h"
  39 #include "AliTPCFast.h"
  40 #include "AliLog.h"
  41
  42 ClassImp(AliTPCFast)
  43   //____________________________________________________________________
  44 AliTPCFast::AliTPCFast(const AliTPCFast &param)
  45               :TObject(param),fParam(0)
  46
  47 {
  48   //
  49   //  copy constructor - dummy
  50   //
  51   fParam = param.fParam;
  52 }
  53 AliTPCFast & AliTPCFast::operator =(const AliTPCFast & param)
  54 {
  55   //
  56   // assignment operator - dummy
  57   //
  58   if (this == &param) return (*this);
  59
  60   fParam=param.fParam;
  61   return (*this);
  62 }
  63
  64 //_____________________________________________________________________________
  65 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
  66 {
  67   //--------------------------------------------------------
  68   // TPC simple cluster generator from hits
  69   // obtained from the TPC Fast Simulator
  70   // The point errors are taken from the parametrization
  71   //--------------------------------------------------------
  72
  73   //-----------------------------------------------------------------
  74   // Origin: Marek Kowalski  IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
  75   //-----------------------------------------------------------------
  76   // Adopted to Marian's cluster data structure by I.Belikov, CERN,
  77   // Jouri.Belikov@cern.ch
  78   //----------------------------------------------------------------
  79
  80   /////////////////////////////////////////////////////////////////////////////
  81   //
  82   //---------------------------------------------------------------------
  83   //   ALICE TPC Cluster Parameters
  84   //--------------------------------------------------------------------
  85
  86
  87
  88   // Cluster width in rphi
  89   const Float_t kACrphi=0.18322;
  90   const Float_t kBCrphi=0.59551e-3;
  91   const Float_t kCCrphi=0.60952e-1;
  92   // Cluster width in z
  93   const Float_t kACz=0.19081;
  94   const Float_t kBCz=0.55938e-3;
  95   const Float_t kCCz=0.30428;
  96
  97
  98   AliLoader* loader = runLoader->GetLoader("TPCLoader");
  99   if (!loader) {
 100     AliError("No TPC loader found");
 101     return;
 102   }
 103   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 104   AliRun* aliRun = runLoader->GetAliRun();
 105   if (!aliRun) {
 106     AliError("Couldn't get AliRun object");
 107     return;
 108   }
 109   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 110   if (!tpc) {
 111     AliError("Couldn't get TPC detector");
 112     return;
 113   }
 114   AliTPCParam* param = tpc->GetParam();
 115   if (!param) {
 116     AliError("No TPC parameters available");
 117     return;
 118   }
 119
 120   //if(fDefaults == 0) SetDefaults();
 121
 122   Float_t sigmaRphi,sigmaZ,clRphi,clZ;
 123   //
 124   TParticle *particle; // pointer to a given particle
 125   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 126   Int_t sector;
 127   Int_t ipart;
 128   Float_t xyz[5];
 129   Float_t pl,pt,tanth,rpad,ratio;
 130   Float_t cph,sph;
 131
 132   //---------------------------------------------------------------
 133   //  Get the access to the tracks
 134   //---------------------------------------------------------------
 135
 136   TTree *tH = loader->TreeH();
 137   if (tH == 0x0){
 138     AliFatal("Can not find TreeH in folder");
 139     return;
 140   }
 141   tpc->SetTreeAddress();
 142
 143   Stat_t ntracks = tH->GetEntries();
 144
 145   //Switch to the output file
 146
 147   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 148
 149   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
 150
 151   runLoader->CdGAFile();
 152   //param->Write(param->GetTitle());
 153
 154   AliTPCClustersArray carray;
 155   carray.Setup(param);
 156   carray.SetClusterType("AliTPCcluster");
 157   carray.MakeTree(loader->TreeR());
 158
 159   Int_t nclusters=0; //cluster counter
 160
 161   //------------------------------------------------------------
 162   // Loop over all sectors (72 sectors for 20 deg
 163   // segmentation for both lower and upper sectors)
 164   // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
 165   // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
 166   //
 167   // First cluster for sector 0 starts at "0"
 168   //------------------------------------------------------------
 169
 170   for(Int_t isec=0;isec<param->GetNSector();isec++){
 171     //MI change
 172     param->AdjustCosSin(isec,cph,sph);
 173
 174     //------------------------------------------------------------
 175     // Loop over tracks
 176     //------------------------------------------------------------
 177
 178     for(Int_t track=0;track<ntracks;track++){
 179       tpc->ResetHits();
 180       tH->GetEvent(track);
 181       //
 182       //  Get number of the TPC hits
 183       //
 184        tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 185
 186       // Loop over hits
 187       //
 188        while(tpcHit){
 189
 190          if (tpcHit->fQ == 0.) {
 191            tpcHit = (AliTPChit*) tpc->NextHit();
 192            continue; //information about track (I.Belikov)
 193          }
 194          sector=tpcHit->fSector; // sector number
 195
 196          if(sector != isec){
 197            tpcHit = (AliTPChit*) tpc->NextHit();
 198            continue;
 199          }
 200          ipart=tpcHit->Track();
 201          particle=aliRun->GetMCApp()->Particle(ipart);
 202          pl=particle->Pz();
 203          pt=particle->Pt();
 204          if(pt < 1.e-9) pt=1.e-9;
 205          tanth=pl/pt;
 206          tanth = TMath::Abs(tanth);
 207          rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
 208          ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
 209
 210          //   space-point resolutions
 211
 212         sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
 213         sigmaZ   =AliTPCcluster::SigmaZ2(rpad,tanth   );
 214
 215         //   cluster widths
 216
 217         clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
 218         clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
 219
 220         // temporary protection
 221
 222         if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
 223         if(sigmaZ < 0.) sigmaZ=0.4e-3;
 224         if(clRphi < 0.) clRphi=2.5e-3;
 225         if(clZ < 0.) clZ=2.5e-5;
 226
 227         //
 228
 229         //
 230         // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
 231         // then the inaccuracy in a X-Y plane is only along Y (pad row)!
 232         //
 233         Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
 234         Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
 235         xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi));   // y
 236           Float_t alpha=(isec < param->GetNInnerSector()) ?
 237           param->GetInnerAngle() : param->GetOuterAngle();
 238           Float_t ymax=xprim*TMath::Tan(0.5*alpha);
 239           if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
 240         xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
 241           if (TMath::Abs(xyz[1])>param->GetZLength(isec)) xyz[1]=tpcHit->Z();
 242         xyz[2]=sigmaRphi;                                     // fSigmaY2
 243         xyz[3]=sigmaZ;                                        // fSigmaZ2
 244         xyz[4]=tpcHit->fQ;                                    // q
 245
 246         AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
 247         if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
 248
 249         Int_t tracks[3]={tpcHit->Track(), -1, -1};
 250         AliTPCcluster cluster(tracks,xyz);
 251
 252         clrow->InsertCluster(&cluster); nclusters++;
 253
 254         tpcHit = (AliTPChit*)tpc->NextHit();
 255
 256
 257       } // end of loop over hits
 258
 259     }   // end of loop over tracks
 260
 261     Int_t nrows=param->GetNRow(isec);
 262     for (Int_t irow=0; irow<nrows; irow++) {
 263         AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
 264         if (!clrow) continue;
 265         carray.StoreRow(isec,irow);
 266         carray.ClearRow(isec,irow);
 267     }
 268
 269   } // end of loop over sectors
 270
 271   //  cerr<<"Number of made clusters : "<<nclusters<<"                        \n";
 272   loader->WriteRecPoints("OVERWRITE");
 273
 274
 275 } // end of function
 276
 277
 278
 279 //_________________________________________________________________
 280
 281
 282
 283 void  AliTPCFast::Hits2ExactClusters(AliRunLoader* runLoader) const{
 284
 285
 286
 287   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 288   if (!loader) {
 289     AliError("No TPC loader found");
 290     return;
 291   }
 292   AliTPC* tpc = (AliTPC*) gAlice->GetDetector("TPC");
 293   if (!tpc) {
 294     AliError("Couldn't get TPC detector");
 295     return;
 296   }
 297   AliTPCParam* param = tpc->GetParam();
 298   if (!param) {
 299     AliError("No TPC parameters available");
 300     return;
 301   }
 302   //---------------------------------------------------------------
 303   //  Get the access to the tracks
 304   //---------------------------------------------------------------
 305
 306   TTree *tH = loader->TreeH();
 307   if (tH == 0x0) { runLoader->LoadHits("TPC"); loader->TreeH();}
 308
 309   tpc->SetTreeAddress();
 310
 311
 312   //Switch to the output file
 313
 314   if (loader->TreeR() == 0x0) loader->MakeTree("R");
 315
 316   AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
 317
 318   runLoader->CdGAFile();
 319   //param->Write(param->GetTitle());
 320
 321   AliTPCClustersArray carray;
 322   carray.Setup(param);
 323   carray.SetClusterType("AliTPCclusterMI");
 324   carray.MakeTree(loader->TreeR());
 325
 326   //------------------------------------------------------------
 327   // Loop over all sectors (72 sectors for 20 deg)
 328   //------------------------------------------------------------
 329
 330   for(Int_t isec=0;isec<param->GetNSector();isec++){
 331     Hits2ExactClustersSector(runLoader, &carray, isec);
 332   } // end of loop over sectors
 333   loader->WriteRecPoints("OVERWRITE");
 334 }
 335
 336
 337
 338
 339 //_________________________________________________________________
 340 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
 341                                           AliTPCClustersArray* clustersArray,
 342                                           Int_t isec) const
 343 {
 344   //--------------------------------------------------------
 345   //calculate exact cross point of track and given pad row
 346   //resulting values are expressed in "local" coordinata
 347   //the sigmaY2 and sigmaZ2 of cluster are the shape of cluster parameters
 348   //                        - thay are used later on for error parameterization
 349   //--------------------------------------------------------
 350
 351   //-----------------------------------------------------------------
 352   // Origin: Marian Ivanov  GSI Darmstadt, m.ivanov@gsi.de
 353   //-----------------------------------------------------------------
 354   //
 355   if (clustersArray==0){
 356     return;
 357   }
 358   AliLoader* loader = runLoader->GetLoader("TPCLoader");
 359   if (!loader) {
 360     AliError("No TPC loader found");
 361     return;
 362   }
 363   if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
 364   AliRun* aliRun = runLoader->GetAliRun();
 365   if (!aliRun) {
 366     AliError("Couldn't get AliRun object");
 367     return;
 368   }
 369   AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
 370   if (!tpc) {
 371     AliError("Couldn't get TPC detector");
 372     return;
 373   }
 374   AliTPCParam* param = tpc->GetParam();
 375   if (!param) {
 376     AliError("No TPC parameters available");
 377     return;
 378   }
 379   //
 380   //
 381   AliTPChit *tpcHit; // pointer to a sigle TPC hit
 382   //  Int_t sector,nhits;
 383   Int_t ipart;
 384   const Int_t kcmaxhits=30000;
 385   TVector * xxxx = new TVector(kcmaxhits*4);
 386   TVector & xxx = *xxxx;
 387   Int_t maxhits = kcmaxhits;
 388   //construct array for each padrow
 389   for (Int_t i=0; i<param->GetNRow(isec);i++)
 390     clustersArray->CreateRow(isec,i);
 391
 392   //---------------------------------------------------------------
 393   //  Get the access to the tracks
 394   //---------------------------------------------------------------
 395
 396   TTree *tH = loader->TreeH();
 397   if (tH == 0x0){
 398     AliFatal("Can not find TreeH in folder");
 399     return;
 400   }
 401   tpc->SetTreeAddress();
 402
 403   Stat_t ntracks = tH->GetEntries();
 404   //MI change
 405   TBranch * branch=0;
 406   if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
 407   else branch = tH->GetBranch("TPC");
 408
 409   //------------------------------------------------------------
 410   // Loop over tracks
 411   //------------------------------------------------------------
 412
 413   for(Int_t track=0;track<ntracks;track++){
 414     Bool_t isInSector=kTRUE;
 415     tpc->ResetHits();
 416     isInSector = tpc->TrackInVolume(isec,track);
 417     if (!isInSector) continue;
 418     //MI change
 419     branch->GetEntry(track); // get next track
 420     //
 421     // Loop over hits
 422     //
 423     Int_t currentIndex=0;
 424     Int_t lastrow=-1;  //last writen row
 425
 426     //M.I. changes
 427
 428     tpcHit = (AliTPChit*)tpc->FirstHit(-1);
 429     while(tpcHit){
 430
 431       Int_t sector=tpcHit->fSector; // sector number
 432       if(sector != isec){
 433         tpcHit = (AliTPChit*) tpc->NextHit();
 434         continue;
 435       }
 436
 437       ipart=tpcHit->Track();
 438
 439       //find row number
 440
 441       Float_t  x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
 442       Int_t    index[3]={1,isec,0};
 443       Int_t    currentrow = param->GetPadRow(x,index);
 444
 445       if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
 446       if (lastrow<0) lastrow=currentrow;
 447       //
 448       //
 449       if (currentrow!=lastrow){
 450         if (currentIndex>2){
 451           Float_t sumx=0;
 452           Float_t sumx2=0;
 453           Float_t sumx3=0;
 454           Float_t sumx4=0;
 455           Float_t sumy=0;
 456           Float_t sumxy=0;
 457           Float_t sumx2y=0;
 458           Float_t sumz=0;
 459           Float_t sumxz=0;
 460           Float_t sumx2z=0;
 461           Float_t sumq=0;
 462           for (Int_t cindex=0;cindex<currentIndex;cindex++){
 463             Float_t x1,x2,x3,x4;
 464             x1=x2=x3=x4=xxx(cindex*4);
 465             x2*=x1;
 466             x3*=x2;
 467             x4*=x3;
 468             sumx+=x1;
 469             sumx2+=x2;
 470             sumx3+=x3;
 471             sumx4+=x4;
 472             sumy+=xxx(cindex*4+1);
 473             sumxy+=xxx(cindex*4+1)*x1;
 474             sumx2y+=xxx(cindex*4+1)*x2;
 475             sumz+=xxx(cindex*4+2);
 476             sumxz+=xxx(cindex*4+2)*x1;
 477             sumx2z+=xxx(cindex*4+2)*x2;
 478             sumq+=xxx(cindex*4+3);
 479           }
 480           Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
 481             sumx2*(sumx*sumx3-sumx2*sumx2);
 482
 483           Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
 484             sumx2*(sumxy*sumx3-sumx2y*sumx2);
 485           Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
 486             sumx2*(sumxz*sumx3-sumx2z*sumx2);
 487
 488           Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
 489             sumx2*(sumx*sumx2y-sumx2*sumxy);
 490           Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
 491             sumx2*(sumx*sumx2z-sumx2*sumxz);
 492
 493           if (TMath::Abs(det)<0.00000000000000001){
 494              tpcHit = (AliTPChit*)tpc->NextHit();
 495             continue;
 496           }
 497
 498           Float_t y=detay/det;
 499           Float_t z=detaz/det;
 500           Float_t by=detby/det; //y angle
 501           Float_t bz=detbz/det; //z angle
 502           sumy/=Float_t(currentIndex);
 503           sumz/=Float_t(currentIndex);
 504
 505
 506           //
 507           //
 508           Float_t sign = (tpcHit->Z()>0)? 1.:-1.;
 509           y = (y+0.5)*param->GetPadPitchWidth(isec);
 510           z = z*param->GetZWidth();
 511           //
 512           // y expected shape
 513           Double_t sigmay2 = z*fParam->GetDiffL();
 514           sigmay2 *= sigmay2;
 515           Float_t angulary = by*param->GetPadPitchLength(isec,lastrow);
 516           angulary*=angulary;
 517           angulary/=12;
 518           sigmay2 +=angulary+0.25*param->GetPadPitchWidth(isec)*param->GetPadPitchWidth(isec);
 519           //
 520           // z expected shape
 521           Double_t sigmaz2 = z*fParam->GetDiffT();
 522           sigmaz2 *= sigmaz2;
 523           Float_t angularz = bz*param->GetPadPitchLength(isec,lastrow);
 524           angularz*=angularz;
 525           angularz/=12;
 526           sigmaz2 +=angularz+0.25;
 527           //
 528           sigmaz2 = TMath::Min(sigmaz2,1.);
 529           sigmay2 = TMath::Min(sigmay2,1.);
 530           //
 531           //
 532           z = sign*(param->GetZLength(isec) - z);
 533           if (TMath::Abs(z)< param->GetZLength(isec)-1){
 534             AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
 535             if (row!=0) {
 536               AliTPCclusterMI* cl = new((AliTPCclusterMI*)row->Append()) AliTPCclusterMI ;
 537               //          AliTPCclusterMI cl;
 538               cl->SetZ(z);
 539               cl->SetY(y);
 540               cl->SetQ(sumq);
 541               if (TMath::Abs(sumq)<1000){
 542                 cl->SetMax(Short_t(sumq));
 543               }
 544               else{
 545                 cl->SetMax(0);
 546               }
 547               cl->SetSigmaZ2(sigmaz2);
 548               cl->SetSigmaY2(sigmay2);
 549               cl->SetLabel(ipart,0);
 550               cl->SetLabel(-1,1);
 551               cl->SetLabel(-1,2);
 552               cl->SetType(0);
 553             }
 554           }
 555         } //end of calculating cluster for given row
 556         currentIndex=0;
 557         lastrow=currentrow;
 558       }  //  end of crossing rows
 559       //
 560       if ( currentIndex>=maxhits){
 561         maxhits+=kcmaxhits;
 562         xxx.ResizeTo(4*maxhits);
 563       }
 564       xxx(currentIndex*4)=x[0];
 565       xxx(currentIndex*4+1)=x[1];
 566       xxx(currentIndex*4+2)=x[2];
 567       xxx(currentIndex*4+3)=tpcHit->fQ;
 568       currentIndex++;
 569       tpcHit = (AliTPChit*)tpc->NextHit();
 570     } // end of loop over hits
 571   }   // end of loop over tracks
 572   //write padrows to tree
 573   for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
 574     clustersArray->StoreRow(isec,ii);
 575     clustersArray->ClearRow(isec,ii);
 576   }
 577   xxxx->Delete();
 578
 579 }
 580
 581