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