TPC/AliBarrelRec_TPCparam.C

   1 /****************************************************************************
   2  * This macro is supposed to do reconstruction in the ITS via Kalman        *
   3  * tracker V2. The ITStracker is feeded with parametrized TPC tracks        *
   4  *                                                                          *
   5  * It does the following steps:                                             *
   6  *             1) TPC tracking parameterization                             *
   7  *             2) ITS cluster finding V2 (via fast points !)                *
   8  *             3) ITS track finding V2                                      *
   9  *             4) Create a reference file with simulation info (p,PDG...)   *
  10  *                                                                          *
  11  * (Origin: A.Dainese, Padova, andrea.dainese@pd,infn.it                    *
  12  *  from AliBarrelReconstruction.C I.Belikov, CERN, Jouri.Belikov@cern.ch)  *
  13  ****************************************************************************/
  14
  15 #ifndef __CINT__
  16   #include "alles.h"
  17   #include "AliMagF.h"
  18   #include "AliITS.h"
  19   #include "AliITSgeom.h"
  20   #include "AliITSRecPoint.h"
  21   #include "AliITSclusterV2.h"
  22   #include "AliITSsimulationFastPoints.h"
  23   #include "AliITStrackerV2.h"
  24   #include "AliTPCtrackerParam.h"
  25 #endif
  26
  27 typedef struct {
  28   Int_t lab;
  29   Int_t pdg;
  30   Int_t mumpdg;
  31   Float_t Vx,Vy,Vz;
  32   Float_t Px,Py,Pz;
  33 } RECTRACK;
  34
  35
  36 Int_t TPCParamTracks(const Char_t *galice,const Char_t *outname,const Int_t coll,const Double_t Bfield,Int_t n);
  37 Int_t ITSFindClusters(const Char_t *inname,const Char_t *outname,Int_t n);
  38 Int_t ITSFindTracks(const Char_t *galice,const Char_t *inname,const Char_t *inname2,const Char_t *outname,Int_t n);
  39 Int_t ITSMakeRefFile(const Char_t *galice, const Char_t *inname, const Char_t *outname, Int_t n);
  40
  41 Int_t AliBarrelRec_TPCparam(Int_t n=1) {
  42   const Char_t *TPCtrkNameS="AliTPCtracksParam.root";
  43   const Char_t *galiceName="galice.root";
  44   const Char_t *ITSclsName="AliITSclustersV2.root";
  45   const Char_t *ITStrkName="AliITStracksV2.root";
  46   const Char_t *ITSrefName="ITStracksRefFile.root";
  47
  48   // set here the code for the type of collision (needed for TPC tracking
  49   // parameterization). available collisions:
  50   //
  51   // coll = 0 ->   PbPb6000 (HIJING with b<2fm)
  52   const Int_t    collcode = 0;
  53   // set here the value of the magnetic field
  54   const Double_t BfieldValue = 0.4;
  55
  56
  57
  58   AliKalmanTrack::SetConvConst(100/0.299792458/BfieldValue);
  59
  60
  61   // ********** Build TPC tracks with parameterization *********** //
  62   if (TPCParamTracks(galiceName,TPCtrkNameS,collcode,BfieldValue,n)) {
  63     cerr<<"Failed to get TPC hits !\n";
  64     return 1;
  65   }
  66
  67
  68   // ********** Find ITS clusters *********** //
  69   if (ITSFindClusters(galiceName,ITSclsName,n)) {
  70     cerr<<"Failed to get ITS clusters !\n";
  71     return 1;
  72   }
  73
  74
  75   // ********* Find ITS tracks *********** //
  76   if (ITSFindTracks(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,n)) {
  77     cerr<<"Failed to get ITS tracks !\n";
  78     return 1;
  79   }
  80
  81
  82   // ********* Make ITS tracks reference file *********** //
  83   if (ITSMakeRefFile(galiceName,ITStrkName,ITSrefName,n)) {
  84     cerr<<"Failed to get ITS tracks ref file!\n";
  85     return 1;
  86   }
  87
  88
  89   return 0;
  90 }
  91
  92
  93 Int_t TPCParamTracks(const Char_t *galice, const Char_t *outname,
  94                      const Int_t coll,const Double_t Bfield,Int_t n) {
  95
  96   cerr<<"\n*******************************************************************\n";
  97   Int_t rc;
  98
  99   const Char_t *name="TPCParamTracks";
 100   cerr<<'\n'<<name<<"...\n";
 101   gBenchmark->Start(name);
 102
 103   TFile *outfile=TFile::Open(outname,"recreate");
 104   TFile *infile =TFile::Open(galice);
 105
 106   AliTPCtrackerParam tracker(coll,Bfield);
 107   rc = tracker.BuildTPCtracks(infile,outfile,n);
 108
 109   delete gAlice; gAlice=0;
 110
 111   infile->Close();
 112   outfile->Close();
 113
 114   gBenchmark->Stop(name);
 115   gBenchmark->Show(name);
 116
 117   return rc;
 118 }
 119
 120 Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname, Int_t n) {
 121
 122
 123   cerr<<"\n*******************************************************************\n";
 124
 125   Int_t rc=0;
 126   const Char_t *name="ITSFindClusters";
 127   cerr<<'\n'<<name<<"...\n";
 128   gBenchmark->Start(name);
 129   TFile *out=TFile::Open(outname,"recreate");
 130   TFile *in =TFile::Open(inname,"update");
 131
 132
 133   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
 134     cerr<<"Can't get gAlice !\n";
 135     return 1;
 136   }
 137
 138
 139   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 140   if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
 141   AliITSgeom *geom=ITS->GetITSgeom();
 142   out->cd();
 143   geom->Write();
 144
 145   Int_t ev=0;
 146   for (ev = 0; ev<n; ev++){
 147     in->cd();   // !!!! MI directory must point to galice. - othervise problem with Tree -connection
 148     gAlice->GetEvent(ev);
 149     //gAlice->TreeR()->Reset();   //reset reconstructed tree
 150
 151
 152     TTree *pTree=gAlice->TreeR();
 153     if (!pTree){
 154       gAlice->MakeTree("R");
 155       pTree = gAlice->TreeR();
 156     }
 157     TBranch *branch=pTree->GetBranch("ITSRecPoints");
 158     if (!branch) {
 159       //if not reconstructed ITS branch do reconstruction
 160       ITS->MakeBranch("R",0);
 161       //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
 162       for (Int_t i=0;i<3;i++) {
 163         ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
 164       }
 165       Int_t nsignal=25;
 166       Int_t size=-1;
 167       Int_t bgr_ev=Int_t(ev/nsignal);
 168       ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
 169       ///////////////////////////////////////////////////////////////////////
 170       gAlice->GetEvent(ev);   //MI comment  - in HitsToFast... they reset treeR to 0
 171       //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
 172       //time
 173     }
 174
 175
 176
 177     out->cd();
 178     TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
 179     char   cname[100];
 180     sprintf(cname,"TreeC_ITS_%d",ev);
 181
 182     TTree *cTree=new TTree(cname,"ITS clusters");
 183     cTree->Branch("Clusters",&clusters);
 184
 185     pTree=gAlice->TreeR();
 186     if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
 187     branch=pTree->GetBranch("ITSRecPoints");
 188     if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
 189     TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
 190     branch->SetAddress(&points);
 191
 192     TClonesArray &cl=*clusters;
 193     Int_t nclusters=0;
 194     Int_t nentr=(Int_t)pTree->GetEntries();
 195     AliITSgeom *geom=ITS->GetITSgeom();
 196
 197     for (Int_t i=0; i<nentr; i++) {
 198       if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
 199       Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
 200       Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
 201       Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
 202       Double_t yshift = x*rot[0] + y*rot[1];
 203       Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
 204       Int_t ncl=points->GetEntriesFast();
 205       nclusters+=ncl;
 206       Float_t lp[5];
 207       Int_t lab[6];
 208       for (Int_t j=0; j<ncl; j++) {
 209         AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
 210         lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
 211         lp[1]=p->GetZ()+zshift;
 212         lp[2]=p->GetSigmaX2();
 213         lp[3]=p->GetSigmaZ2();
 214         lp[4]=p->GetQ();
 215         lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
 216         lab[3]=ndet;
 217
 218         Int_t label=lab[0];
 219         TParticle *part=(TParticle*)gAlice->Particle(label);
 220         label=-3;
 221         while (part->P() < 0.005) {
 222           Int_t m=part->GetFirstMother();
 223           if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
 224           label=m;
 225           part=(TParticle*)gAlice->Particle(label);
 226         }
 227         if      (lab[1]<0) lab[1]=label;
 228         else if (lab[2]<0) lab[2]=label;
 229         else cerr<<"No empty labels !\n";
 230
 231         new(cl[j]) AliITSclusterV2(lab,lp);
 232       }
 233       cTree->Fill(); clusters->Delete();
 234       points->Delete();
 235     }
 236     cTree->Write();
 237     cerr<<"Number of clusters: "<<nclusters<<endl;
 238     delete cTree; delete clusters; delete points;
 239
 240   }
 241
 242
 243   delete gAlice; gAlice=0;
 244   in->Close();
 245   out->Close();
 246   gBenchmark->Stop(name);
 247   gBenchmark->Show(name);
 248
 249   return rc;
 250 }
 251
 252 Int_t ITSFindTracks(const Char_t *galice, const Char_t * inname,
 253                     const Char_t *inname2, const Char_t *outname,
 254                     Int_t n) {
 255
 256
 257   cerr<<"\n*******************************************************************\n";
 258
 259   Int_t rc=0;
 260   const Char_t *name="ITSFindTracks";
 261   cerr<<'\n'<<name<<"...\n";
 262   gBenchmark->Start(name);
 263
 264
 265   TFile *out=TFile::Open(outname,"recreate");
 266   TFile *in =TFile::Open(inname);
 267   TFile *in2 =TFile::Open(inname2);
 268
 269   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
 270   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
 271
 272
 273   for (Int_t ev=0; ev<n; ev++){
 274     AliITStrackerV2 tracker(geom,ev);
 275
 276     TArrayF o;
 277     o.Set(3);
 278     TFile* vtxFile = new TFile(galice);
 279     vtxFile->cd();
 280     AliHeader* header = 0;
 281     TTree* treeE = (TTree*)gDirectory->Get("TE");
 282     treeE->SetBranchAddress("Header",&header);
 283     treeE->GetEntry(ev);
 284     AliGenEventHeader* genHeader = header->GenEventHeader();
 285     if(genHeader) {
 286       // get primary vertex position
 287       genHeader->PrimaryVertex(o);
 288       vtxFile->Close();
 289       delete header;
 290     // set position of primary vertex
 291       Double_t vtx[3];
 292       vtx[0]=o[0]; vtx[1]=o[1]; vtx[2]=o[2];
 293       cerr<<"+++\n+++ Reading primary vertex position from galice.root\n+++\n+++ Setting primary vertex z = "<<vtx[2]<<" cm for ITS tracking\n+++\n";
 294       tracker.SetVertex(vtx);
 295     }else {
 296       cerr<<"+++\n+++ Event header not found in galice.root:\n+++ Primary vertex in (0,0,0) [default]\n+++\n";
 297     }
 298
 299     // setup vertex constraint in the two tracking passes
 300     Int_t flags[2];
 301     flags[0]=0;
 302     tracker.SetupFirstPass(flags);
 303     flags[0]=-1;
 304     tracker.SetupSecondPass(flags);
 305
 306     rc=tracker.Clusters2Tracks(in,out);
 307
 308   }
 309
 310   delete gAlice;  gAlice=0;
 311   in->Close();
 312   in2->Close();
 313   out->Close();
 314
 315   gBenchmark->Stop(name);
 316   gBenchmark->Show(name);
 317
 318   return rc;
 319 }
 320
 321
 322 Int_t ITSMakeRefFile(const Char_t *galice, const Char_t *inname,
 323                      const Char_t *outname, Int_t n) {
 324
 325
 326   cerr<<"\n*******************************************************************\n";
 327
 328   Int_t rc=0;
 329   const Char_t *name="ITSMakeRefFile";
 330   cerr<<'\n'<<name<<"...\n";
 331   gBenchmark->Start(name);
 332
 333
 334   TFile *out = TFile::Open(outname,"recreate");
 335   TFile *trk = TFile::Open(inname);
 336   TFile *kin = TFile::Open(galice);
 337
 338
 339   // Get gAlice object from file
 340   if(!(gAlice=(AliRun*)kin->Get("gAlice"))) {
 341     cerr<<"gAlice has not been found on galice.root !\n";
 342     return 1;
 343   }
 344
 345
 346
 347   Int_t label;
 348   TParticle *Part;
 349   TParticle *Mum;
 350   static RECTRACK rectrk;
 351
 352
 353   for(Int_t event=0; event<n; event++){
 354
 355     AliITStrackV2 *itstrack=0;
 356
 357     Int_t nparticles=gAlice->GetEvent(event);
 358
 359     trk->cd();
 360
 361     // Tree with ITS tracks
 362     char tname[100];
 363     sprintf(tname,"TreeT_ITS_%d",event);
 364
 365     TTree *tracktree=(TTree*)trk->Get(tname);
 366     TBranch *tbranch=tracktree->GetBranch("tracks");
 367     Int_t nentr=(Int_t)tracktree->GetEntries();
 368
 369     // Tree for true track parameters
 370     char ttname[100];
 371     sprintf(ttname,"Tree_Ref_%d",event);
 372     TTree *reftree = new TTree(ttname,"Tree with true track params");
 373     reftree->Branch("rectracks",&rectrk,"lab/I:pdg:Vx/F:Vy:Vz:Px:Py:Pz");
 374
 375     for (Int_t i=0; i<nentr; i++) {
 376       itstrack=new AliITStrackV2;
 377       tbranch->SetAddress(&itstrack);
 378       tracktree->GetEvent(i);
 379       label = TMath::Abs(itstrack->GetLabel());
 380
 381       Part = (TParticle*)gAlice->Particle(label);
 382       rectrk.lab=label;
 383       rectrk.pdg=Part->GetPdgCode();
 384       if(Part->GetFirstMother()>=0) {
 385         Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother());
 386         rectrk.mumpdg=Mum->GetPdgCode();
 387       } else {
 388         rectrk.mumpdg=-1;
 389       }
 390       rectrk.Vx=Part->Vx();
 391       rectrk.Vy=Part->Vy();
 392       rectrk.Vz=Part->Vz();
 393       rectrk.Px=Part->Px();
 394       rectrk.Py=Part->Py();
 395       rectrk.Pz=Part->Pz();
 396
 397       reftree->Fill();
 398     }
 399
 400     out->cd();
 401     reftree->Write();
 402
 403   }
 404
 405   trk->Close();
 406   kin->Close();
 407   out->Close();
 408
 409   gBenchmark->Stop(name);
 410   gBenchmark->Show(name);
 411
 412
 413   return rc;
 414
 415 }
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427