1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////////////////
20 // fast TPC cluster simulation //
22 ///////////////////////////////////////////////////////////////////////////////
24 #include <TParticle.h>
28 #include "AliRunLoader.h"
32 #include "AliTPCParam.h"
33 #include "AliTPCClustersArray.h"
34 #include "AliTPCClustersRow.h"
35 #include "AliTPCcluster.h"
36 #include "AliComplexCluster.h"
37 #include "AliTPCFast.h"
41 //____________________________________________________________________
42 AliTPCFast::AliTPCFast(const AliTPCFast ¶m)
43 :TObject(param),fParam(0)
47 // copy constructor - dummy
49 fParam = param.fParam;
51 AliTPCFast & AliTPCFast::operator =(const AliTPCFast & param)
54 // assignment operator - dummy
60 //_____________________________________________________________________________
61 void AliTPCFast::Hits2Clusters(AliRunLoader* runLoader) const
63 //--------------------------------------------------------
64 // TPC simple cluster generator from hits
65 // obtained from the TPC Fast Simulator
66 // The point errors are taken from the parametrization
67 //--------------------------------------------------------
69 //-----------------------------------------------------------------
70 // Origin: Marek Kowalski IFJ, Krakow, Marek.Kowalski@ifj.edu.pl
71 //-----------------------------------------------------------------
72 // Adopted to Marian's cluster data structure by I.Belikov, CERN,
73 // Jouri.Belikov@cern.ch
74 //----------------------------------------------------------------
76 /////////////////////////////////////////////////////////////////////////////
78 //---------------------------------------------------------------------
79 // ALICE TPC Cluster Parameters
80 //--------------------------------------------------------------------
84 // Cluster width in rphi
85 const Float_t kACrphi=0.18322;
86 const Float_t kBCrphi=0.59551e-3;
87 const Float_t kCCrphi=0.60952e-1;
89 const Float_t kACz=0.19081;
90 const Float_t kBCz=0.55938e-3;
91 const Float_t kCCz=0.30428;
94 AliLoader* loader = runLoader->GetLoader("TPCLoader");
96 AliError("No TPC loader found");
99 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
100 AliRun* aliRun = runLoader->GetAliRun();
102 AliError("Couldn't get AliRun object");
105 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
107 AliError("Couldn't get TPC detector");
110 AliTPCParam* param = tpc->GetParam();
112 AliError("No TPC parameters available");
116 //if(fDefaults == 0) SetDefaults();
118 Float_t sigmaRphi,sigmaZ,clRphi,clZ;
120 TParticle *particle; // pointer to a given particle
121 AliTPChit *tpcHit; // pointer to a sigle TPC hit
125 Float_t pl,pt,tanth,rpad,ratio;
128 //---------------------------------------------------------------
129 // Get the access to the tracks
130 //---------------------------------------------------------------
132 TTree *tH = loader->TreeH();
134 AliFatal("Can not find TreeH in folder");
136 tpc->SetTreeAddress();
138 Stat_t ntracks = tH->GetEntries();
140 //Switch to the output file
142 if (loader->TreeR() == 0x0) loader->MakeTree("R");
144 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
146 runLoader->CdGAFile();
147 //param->Write(param->GetTitle());
149 AliTPCClustersArray carray;
151 carray.SetClusterType("AliTPCcluster");
152 carray.MakeTree(loader->TreeR());
154 Int_t nclusters=0; //cluster counter
156 //------------------------------------------------------------
157 // Loop over all sectors (72 sectors for 20 deg
158 // segmentation for both lower and upper sectors)
159 // Sectors 0-35 are lower sectors, 0-17 z>0, 17-35 z<0
160 // Sectors 36-71 are upper sectors, 36-53 z>0, 54-71 z<0
162 // First cluster for sector 0 starts at "0"
163 //------------------------------------------------------------
165 for(Int_t isec=0;isec<param->GetNSector();isec++){
167 param->AdjustCosSin(isec,cph,sph);
169 //------------------------------------------------------------
171 //------------------------------------------------------------
173 for(Int_t track=0;track<ntracks;track++){
177 // Get number of the TPC hits
179 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
185 if (tpcHit->fQ == 0.) {
186 tpcHit = (AliTPChit*) tpc->NextHit();
187 continue; //information about track (I.Belikov)
189 sector=tpcHit->fSector; // sector number
192 tpcHit = (AliTPChit*) tpc->NextHit();
195 ipart=tpcHit->Track();
196 particle=aliRun->GetMCApp()->Particle(ipart);
199 if(pt < 1.e-9) pt=1.e-9;
201 tanth = TMath::Abs(tanth);
202 rpad=TMath::Sqrt(tpcHit->X()*tpcHit->X() + tpcHit->Y()*tpcHit->Y());
203 ratio=0.001*rpad/pt; // pt must be in MeV/c - historical reason
205 // space-point resolutions
207 sigmaRphi=AliTPCcluster::SigmaY2(rpad,tanth,pt);
208 sigmaZ =AliTPCcluster::SigmaZ2(rpad,tanth );
212 clRphi=kACrphi-kBCrphi*rpad*tanth+kCCrphi*ratio*ratio;
213 clZ=kACz-kBCz*rpad*tanth+kCCz*tanth*tanth;
215 // temporary protection
217 if(sigmaRphi < 0.) sigmaRphi=0.4e-3;
218 if(sigmaZ < 0.) sigmaZ=0.4e-3;
219 if(clRphi < 0.) clRphi=2.5e-3;
220 if(clZ < 0.) clZ=2.5e-5;
225 // smearing --> rotate to the 1 (13) or to the 25 (49) sector,
226 // then the inaccuracy in a X-Y plane is only along Y (pad row)!
228 Float_t xprim= tpcHit->X()*cph + tpcHit->Y()*sph;
229 Float_t yprim=-tpcHit->X()*sph + tpcHit->Y()*cph;
230 xyz[0]=gRandom->Gaus(yprim,TMath::Sqrt(sigmaRphi)); // y
231 Float_t alpha=(isec < param->GetNInnerSector()) ?
232 param->GetInnerAngle() : param->GetOuterAngle();
233 Float_t ymax=xprim*TMath::Tan(0.5*alpha);
234 if (TMath::Abs(xyz[0])>ymax) xyz[0]=yprim;
235 xyz[1]=gRandom->Gaus(tpcHit->Z(),TMath::Sqrt(sigmaZ)); // z
236 if (TMath::Abs(xyz[1])>param->GetZLength(isec)) xyz[1]=tpcHit->Z();
237 xyz[2]=sigmaRphi; // fSigmaY2
238 xyz[3]=sigmaZ; // fSigmaZ2
239 xyz[4]=tpcHit->fQ; // q
241 AliTPCClustersRow *clrow=carray.GetRow(sector,tpcHit->fPadRow);
242 if (!clrow) clrow=carray.CreateRow(sector,tpcHit->fPadRow);
244 Int_t tracks[3]={tpcHit->Track(), -1, -1};
245 AliTPCcluster cluster(tracks,xyz);
247 clrow->InsertCluster(&cluster); nclusters++;
249 tpcHit = (AliTPChit*)tpc->NextHit();
252 } // end of loop over hits
254 } // end of loop over tracks
256 Int_t nrows=param->GetNRow(isec);
257 for (Int_t irow=0; irow<nrows; irow++) {
258 AliTPCClustersRow *clrow=carray.GetRow(isec,irow);
259 if (!clrow) continue;
260 carray.StoreRow(isec,irow);
261 carray.ClearRow(isec,irow);
264 } // end of loop over sectors
266 // cerr<<"Number of made clusters : "<<nclusters<<" \n";
267 loader->WriteRecPoints("OVERWRITE");
274 //_________________________________________________________________
278 void AliTPCFast::Hits2ExactClusters(AliRunLoader* runLoader) const{
282 AliLoader* loader = runLoader->GetLoader("TPCLoader");
284 AliError("No TPC loader found");
287 AliTPC* tpc = (AliTPC*) gAlice->GetDetector("TPC");
289 AliError("Couldn't get TPC detector");
292 AliTPCParam* param = tpc->GetParam();
294 AliError("No TPC parameters available");
297 //---------------------------------------------------------------
298 // Get the access to the tracks
299 //---------------------------------------------------------------
301 TTree *tH = loader->TreeH();
302 if (tH == 0x0) { runLoader->LoadHits("TPC"); loader->TreeH();}
304 tpc->SetTreeAddress();
307 //Switch to the output file
309 if (loader->TreeR() == 0x0) loader->MakeTree("R");
311 AliDebug(1,Form("param->GetTitle() = %s",param->GetTitle()));
313 runLoader->CdGAFile();
314 //param->Write(param->GetTitle());
316 AliTPCClustersArray carray;
318 carray.SetClusterType("AliTPCclusterMI");
319 carray.MakeTree(loader->TreeR());
321 //------------------------------------------------------------
322 // Loop over all sectors (72 sectors for 20 deg)
323 //------------------------------------------------------------
325 for(Int_t isec=0;isec<param->GetNSector();isec++){
326 Hits2ExactClustersSector(runLoader, &carray, isec);
327 } // end of loop over sectors
328 loader->WriteRecPoints("OVERWRITE");
334 //_________________________________________________________________
335 void AliTPCFast::Hits2ExactClustersSector(AliRunLoader* runLoader,
336 AliTPCClustersArray* clustersArray,
339 //--------------------------------------------------------
340 //calculate exact cross point of track and given pad row
341 //resulting values are expressed in "local" coordinata
342 //the sigmaY2 and sigmaZ2 of cluster are the shape of cluster parameters
343 // - thay are used later on for error parameterization
344 //--------------------------------------------------------
346 //-----------------------------------------------------------------
347 // Origin: Marian Ivanov GSI Darmstadt, m.ivanov@gsi.de
348 //-----------------------------------------------------------------
350 if (clustersArray==0){
353 AliLoader* loader = runLoader->GetLoader("TPCLoader");
355 AliError("No TPC loader found");
358 if (!runLoader->GetAliRun()) runLoader->LoadgAlice();
359 AliRun* aliRun = runLoader->GetAliRun();
361 AliError("Couldn't get AliRun object");
364 AliTPC* tpc = (AliTPC*) aliRun->GetDetector("TPC");
366 AliError("Couldn't get TPC detector");
369 AliTPCParam* param = tpc->GetParam();
371 AliError("No TPC parameters available");
376 AliTPChit *tpcHit; // pointer to a sigle TPC hit
377 // Int_t sector,nhits;
379 const Int_t kcmaxhits=30000;
380 TVector * xxxx = new TVector(kcmaxhits*4);
381 TVector & xxx = *xxxx;
382 Int_t maxhits = kcmaxhits;
383 //construct array for each padrow
384 for (Int_t i=0; i<param->GetNRow(isec);i++)
385 clustersArray->CreateRow(isec,i);
387 //---------------------------------------------------------------
388 // Get the access to the tracks
389 //---------------------------------------------------------------
391 TTree *tH = loader->TreeH();
393 AliFatal("Can not find TreeH in folder");
395 tpc->SetTreeAddress();
397 Stat_t ntracks = tH->GetEntries();
400 if (tpc->GetHitType()>1) branch = tH->GetBranch("TPC2");
401 else branch = tH->GetBranch("TPC");
403 //------------------------------------------------------------
405 //------------------------------------------------------------
407 for(Int_t track=0;track<ntracks;track++){
408 Bool_t isInSector=kTRUE;
410 isInSector = tpc->TrackInVolume(isec,track);
411 if (!isInSector) continue;
413 branch->GetEntry(track); // get next track
417 Int_t currentIndex=0;
418 Int_t lastrow=-1; //last writen row
422 tpcHit = (AliTPChit*)tpc->FirstHit(-1);
425 Int_t sector=tpcHit->fSector; // sector number
427 tpcHit = (AliTPChit*) tpc->NextHit();
431 ipart=tpcHit->Track();
435 Float_t x[3]={tpcHit->X(),tpcHit->Y(),tpcHit->Z()};
436 Int_t index[3]={1,isec,0};
437 Int_t currentrow = param->GetPadRow(x,index);
439 if (currentrow<0) {tpcHit = (AliTPChit*)tpc->NextHit(); continue;}
440 if (lastrow<0) lastrow=currentrow;
443 if (currentrow!=lastrow){
456 for (Int_t cindex=0;cindex<currentIndex;cindex++){
458 x1=x2=x3=x4=xxx(cindex*4);
466 sumy+=xxx(cindex*4+1);
467 sumxy+=xxx(cindex*4+1)*x1;
468 sumx2y+=xxx(cindex*4+1)*x2;
469 sumz+=xxx(cindex*4+2);
470 sumxz+=xxx(cindex*4+2)*x1;
471 sumx2z+=xxx(cindex*4+2)*x2;
472 sumq+=xxx(cindex*4+3);
474 Float_t det=currentIndex*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumx*sumx4-sumx2*sumx3)+
475 sumx2*(sumx*sumx3-sumx2*sumx2);
477 Float_t detay=sumy*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxy*sumx4-sumx2y*sumx3)+
478 sumx2*(sumxy*sumx3-sumx2y*sumx2);
479 Float_t detaz=sumz*(sumx2*sumx4-sumx3*sumx3)-sumx*(sumxz*sumx4-sumx2z*sumx3)+
480 sumx2*(sumxz*sumx3-sumx2z*sumx2);
482 Float_t detby=currentIndex*(sumxy*sumx4-sumx2y*sumx3)-sumy*(sumx*sumx4-sumx2*sumx3)+
483 sumx2*(sumx*sumx2y-sumx2*sumxy);
484 Float_t detbz=currentIndex*(sumxz*sumx4-sumx2z*sumx3)-sumz*(sumx*sumx4-sumx2*sumx3)+
485 sumx2*(sumx*sumx2z-sumx2*sumxz);
487 if (TMath::Abs(det)<0.00000000000000001){
488 tpcHit = (AliTPChit*)tpc->NextHit();
494 Float_t by=detby/det; //y angle
495 Float_t bz=detbz/det; //z angle
496 sumy/=Float_t(currentIndex);
497 sumz/=Float_t(currentIndex);
502 Float_t sign = (tpcHit->Z()>0)? 1.:-1.;
503 y = (y+0.5)*param->GetPadPitchWidth(isec);
504 z = z*param->GetZWidth();
507 Double_t sigmay2 = z*fParam->GetDiffL();
509 Float_t angulary = by*param->GetPadPitchLength(isec,lastrow);
512 sigmay2 +=angulary+0.25*param->GetPadPitchWidth(isec)*param->GetPadPitchWidth(isec);
515 Double_t sigmaz2 = z*fParam->GetDiffT();
517 Float_t angularz = bz*param->GetPadPitchLength(isec,lastrow);
520 sigmaz2 +=angularz+0.25;
522 sigmaz2 = TMath::Min(sigmaz2,1.);
523 sigmay2 = TMath::Min(sigmay2,1.);
526 z = sign*(param->GetZLength(isec) - z);
527 if (TMath::Abs(z)< param->GetZLength(isec)-1){
528 AliTPCClustersRow * row = (clustersArray->GetRow(isec,lastrow));
530 AliTPCclusterMI* cl = new((AliTPCclusterMI*)row->Append()) AliTPCclusterMI ;
531 // AliTPCclusterMI cl;
535 if (TMath::Abs(sumq)<1000){
536 cl->SetMax(Short_t(sumq));
541 cl->SetSigmaZ2(sigmaz2);
542 cl->SetSigmaY2(sigmay2);
543 cl->SetLabel(ipart,0);
549 } //end of calculating cluster for given row
552 } // end of crossing rows
554 if ( currentIndex>=maxhits){
556 xxx.ResizeTo(4*maxhits);
558 xxx(currentIndex*4)=x[0];
559 xxx(currentIndex*4+1)=x[1];
560 xxx(currentIndex*4+2)=x[2];
561 xxx(currentIndex*4+3)=tpcHit->fQ;
563 tpcHit = (AliTPChit*)tpc->NextHit();
564 } // end of loop over hits
565 } // end of loop over tracks
566 //write padrows to tree
567 for (Int_t ii=0; ii<param->GetNRow(isec);ii++) {
568 clustersArray->StoreRow(isec,ii);
569 clustersArray->ClearRow(isec,ii);