TPC/AliBarrelRec_TPCparam.C

   1 /****************************************************************************
   2  * This macro performs track and vertex reconstruction in TPC and ITS.      *
   3  * The ITS Kalman tracker V2 is feeded "with" parameterized TPC tracks.     *
   4  *                                                                          *
   5  * Reconstruction is performed in the following steps:                      *
   6  *             1) TPC tracking parameterization                             *
   7  *             2) ITS clusters: slow or fast                                *
   8  *             3) Primary vertex reconstruction                             *
   9  *                - read from event header for Pb-Pb events                 *
  10  *                - determined using points in pixels for pp/pA events      *
  11  *             4) ITS track finding V2                                      *
  12  *                - in pp/pA, redetermine the position of primary vertex    *
  13  *                  using the reconstructed tracks                          *
  14  *             5) Create a reference file with simulation info (p,PDG...)   *
  15  *                                                                          *
  16  * If mode='A' all 5 steps are executed                                     *
  17  * If mode='B' only steps 4-5 are executed                                  *
  18  *                                                                          *
  19  *  Origin: A.Dainese, Padova,   andrea.dainese@pd.infn.it                  *
  20  *  (from AliTPCtest.C & AliITStestV2.C by I.Belikov)                       *
  21  ****************************************************************************/
  22
  23 // structure for track references
  24 typedef struct {
  25   Int_t lab;
  26   Int_t pdg;
  27   Int_t mumlab;
  28   Int_t mumpdg;
  29   Float_t Vx,Vy,Vz;
  30   Float_t Px,Py,Pz;
  31 } RECTRACK;
  32
  33 //===== Functions definition =================================================
  34
  35 void CopyVtx(const Char_t *inName,const Char_t *outName);
  36
  37 void ITSFindClustersV2(Char_t SlowOrFast);
  38
  39 void ITSFindTracksV2(Int_t *skipEvt);
  40
  41 void ITSMakeRefFile(Int_t *skipEvt);
  42
  43 void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
  44
  45 void PrimaryVertex(const Char_t *outName,Char_t vtxMode);
  46
  47 void TPCParamTracks(Int_t coll,Double_t Bfield);
  48
  49 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName);
  50
  51 void VtxFromHeader(const Char_t *outName,Bool_t smear);
  52
  53 void VtxFromTracks(const Char_t *outName);
  54
  55 void ZvtxFromSPD(const Char_t *outName);
  56
  57 //=============================================================================
  58
  59 // number of events to be processed
  60 Int_t    gNevents;
  61 // magnetic field
  62 Double_t gBfieldValue;
  63
  64 void AliBarrelRec_TPCparam(Int_t n=-1,Char_t mode='A') {
  65
  66   //---------------------------------------------------------------------
  67   //                    CONFIGURATION
  68   //
  69   // _Magnetic_field_
  70   gBfieldValue = 0.4;
  71   //
  72   // _Type_of_collision_ (needed for TPC tracking parameterization)
  73   // Available choices:   !!! ONLY B = 0.4 TESLA !!!
  74   //    collcode = 0  ->   PbPb6000 (HIJING with b<2fm)
  75   //    collcode = 1  ->   low multiplicity: pp or pA
  76   Int_t collcode = 1;
  77   //
  78   // _ITS_clusters_reconstruction_
  79   // Available choices:  (from AliITStestV2.C)
  80   //    SlowOrFast = 's'    slow points
  81   //    SlowOrFast = 'f'    fast points
  82   Char_t SlowOrFast = 'f';
  83   //
  84   // _Primary_vertex_for_ITS_tracking_
  85   // Available choices:
  86   //    Vtx4Tracking = 'H'   from event Header
  87   //    --- for Pb-Pb ---
  88   //    Vtx4Tracking = 'S'   from event header + Smearing
  89   //                                           (x=15,y=15,z=10) micron
  90   //    --- for pp/pA ---
  91   //    Vtx4Tracking = 'P'   z from pixels, x,y in(0,0)
  92   Char_t Vtx4Tracking = 'P';
  93   // _Primary_vertex_for_analysis_ (AliITSVertex stored in tracks file)
  94   // Available choices:
  95   //    Vtx4Analysis = 'C'   Copy the same used for tracking
  96   //    --- for pp/pA ---
  97   //    Vtx4Analysis = 'T'   x,y,z from Tracks
  98   Char_t Vtx4Analysis = 'T';
  99   //
 100   //                  END CONFIGURATION
 101   //---------------------------------------------------------------------
 102
 103   const Char_t *name=" AliBarrelRec_TPCparam";
 104   printf("\n %s\n",name);
 105   gBenchmark->Start(name);
 106
 107   if(n==-1) { // read number of events to be processed from file
 108     TFile *f = new TFile("galice.root");
 109     gAlice = (AliRun*)f->Get("gAlice");
 110     n = gAlice->GetEventsPerRun();
 111     delete gAlice;
 112     gAlice=0;
 113     f->Close();
 114     delete f;
 115     printf(" All %d events in file will be processed\n",n);
 116   }
 117   gNevents = n;
 118
 119
 120   // conversion constant for kalman tracks
 121   AliKalmanTrack::SetConvConst(100/0.299792458/gBfieldValue);
 122
 123   // Selection of execution mode
 124   switch(mode) {
 125   case 'A':
 126     // Build TPC tracks with parameterization
 127     TPCParamTracks(collcode,gBfieldValue);
 128
 129     // ITS clusters
 130     ITSFindClustersV2(SlowOrFast);
 131
 132     // Vertex for ITS tracking
 133     PrimaryVertex("Vtx4Tracking.root",Vtx4Tracking);
 134
 135     break;
 136
 137   case 'B':
 138     printf("       ---> only tracking in ITS <---\n");
 139
 140     // Update list of events to be skipped
 141     if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat")) return;
 142
 143     break;
 144   }
 145
 146   // Mark events that have to be skipped (if any)
 147   Int_t *skipEvt = new Int_t[gNevents];
 148   for(Int_t i=0; i<gNevents; i++) skipEvt[i] = 0;
 149   if(!gSystem->AccessPathName("evtsToSkip.dat",kFileExists))
 150     MarkEvtsToSkip("evtsToSkip.dat",skipEvt);
 151
 152   // Tracking in ITS
 153   ITSFindTracksV2(skipEvt);
 154
 155   // Vertex for analysis
 156   PrimaryVertex("AliITStracksV2.root",Vtx4Analysis);
 157
 158   // Create ITS tracks reference file
 159   ITSMakeRefFile(skipEvt);
 160   delete [] skipEvt;
 161
 162
 163
 164   gBenchmark->Stop(name);
 165   gBenchmark->Show(name);
 166
 167   return;
 168 }
 169 //-----------------------------------------------------------------------------
 170 void CopyVtx(const Char_t *inName,const Char_t *outName) {
 171
 172   // Open input and output files
 173   TFile *inFile = new TFile(inName);
 174   TFile *outFile = new TFile(outName,"update");
 175
 176   TDirectory *curdir;
 177   Char_t vname[20];
 178
 179
 180   for(Int_t ev=0; ev<gNevents; ev++) {
 181     sprintf(vname,"Vertex_%d",ev);
 182     AliITSVertex *vertex = (AliITSVertex*)inFile->Get(vname);
 183     if(!vertex) continue;
 184     curdir = gDirectory;
 185     outFile->cd();
 186     vertex->Write();
 187     curdir->cd();
 188     vertex = 0;
 189   }
 190
 191   inFile->Close();
 192   outFile->Close();
 193   delete inFile;
 194   delete outFile;
 195
 196   return;
 197 }
 198 //-----------------------------------------------------------------------------
 199 void ITSFindClustersV2(Char_t SlowOrFast) {
 200
 201   printf("\n------------------------------------\n");
 202
 203   const Char_t *name="ITSFindClustersV2";
 204   printf("\n %s\n",name);
 205   gBenchmark->Start(name);
 206
 207   //---  taken from AliITStestV2.C--------------------------------------
 208   //
 209   if (SlowOrFast=='f') {
 210     //cerr<<"Fast AliITSRecPoint(s) !\n";
 211     //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
 212     //AliITSHits2FastRecPoints();
 213   } else {
 214     gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2SDigits.C");
 215     AliITSHits2SDigits();
 216     gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSSDigits2Digits.C");
 217     AliITSSDigits2Digits();
 218     //gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSDigits2RecPoints.C");
 219     //AliITSDigits2RecPoints();
 220   }
 221   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSFindClustersV2.C");
 222   AliITSFindClustersV2(SlowOrFast,gNevents);
 223   //
 224   //--------------------------------------------------------------------
 225
 226
 227   gBenchmark->Stop(name);
 228   gBenchmark->Show(name);
 229
 230    return;
 231 }
 232 //-----------------------------------------------------------------------------
 233 Int_t ITSFindTracksV2(Int_t *skipEvt) {
 234
 235   printf("\n------------------------------------\n");
 236
 237   const Char_t *name="ITSFindTracksV2";
 238   printf("\n %s\n",name);
 239   gBenchmark->Start(name);
 240
 241
 242   TFile *outFile     = new TFile("AliITStracksV2.root","recreate");
 243   TFile *inTPCtrks   = new TFile("AliTPCtracksParam.root");
 244   TFile *inVertex    = new TFile("Vtx4Tracking.root");
 245   TFile *inClusters  = new TFile("AliITSclustersV2.root");
 246
 247   AliITSgeom *geom=(AliITSgeom*)inClusters->Get("AliITSgeom");
 248   if(!geom) { printf("can't get ITS geometry !\n"); return;}
 249
 250   Double_t vtx[3];
 251   Int_t flag1stPass,flag2ndPass;
 252   Char_t vname[20];
 253
 254   // open logfile for done events
 255   FILE *logfile = fopen("itstracking.log","w");
 256
 257   // Instantiate AliITStrackerV2
 258   AliITStrackerV2 tracker(geom);
 259
 260   // loop on events
 261   for(Int_t ev=0; ev<gNevents; ev++){
 262     // write to logfile of begun events
 263     fprintf(logfile,"%d\n",ev);
 264
 265     if(skipEvt[ev]) continue;
 266     printf(" --- Processing event %d ---\n",ev);
 267
 268     // pass event number to the tracker
 269     tracker.SetEventNumber(ev);
 270
 271     // set position of primary vertex
 272     sprintf(vname,"Vertex_%d",ev);
 273     AliITSVertex *vertex = (AliITSVertex*)inVertex->Get(vname);
 274     if(vertex) {
 275       vertex->GetXYZ(vtx);
 276       delete vertex;
 277     } else {
 278       printf(" AliITSVertex not found for event %d\n",ev);
 279       printf(" Using (0,0,0) for ITS tracking\n");
 280       vtx[0] = vtx[1] = vtx[2] = 0.;
 281     }
 282
 283     flag1stPass=1; // vtx constraint
 284     flag2ndPass=0; // no vtx constraint
 285
 286     // no vtx constraint if vertex not found
 287     if(vtx[2]<-999.) {
 288       flag1stPass=0;
 289       vtx[2]=0.;
 290     }
 291
 292     tracker.SetVertex(vtx);
 293
 294     // setup vertex constraint in the two tracking passes
 295     Int_t flags[2];
 296     flags[0]=flag1stPass;
 297     tracker.SetupFirstPass(flags);
 298     flags[0]=flag2ndPass;
 299     tracker.SetupSecondPass(flags);
 300
 301     // find the tracks
 302     tracker.Clusters2Tracks(inTPCtrks,outFile);
 303
 304   } // loop on events
 305
 306   fprintf(logfile,"%d\n",gNevents); //this means all evts are successfully completed
 307   fclose(logfile);
 308
 309   delete geom;
 310
 311   inTPCtrks->Close();
 312   inClusters->Close();
 313   inVertex->Close();
 314   outFile->Close();
 315
 316
 317   gBenchmark->Stop(name);
 318   gBenchmark->Show(name);
 319
 320   return;
 321 }
 322 //-----------------------------------------------------------------------------
 323 void ITSMakeRefFile(Int_t *skipEvt) {
 324
 325   printf("\n------------------------------------\n");
 326
 327   const Char_t *name="ITSMakeRefFile";
 328   printf("\n %s\n",name);
 329   gBenchmark->Start(name);
 330
 331
 332   TFile *out = TFile::Open("ITStracksRefFile.root","recreate");
 333   TFile *trk = TFile::Open("AliITStracksV2.root");
 334   TFile *kin = TFile::Open("galice.root");
 335
 336
 337   // Get gAlice object from file
 338   gAlice=(AliRun*)kin->Get("gAlice");
 339
 340   Int_t label;
 341   TParticle *Part;
 342   TParticle *Mum;
 343   RECTRACK rectrk;
 344
 345
 346   for(Int_t ev=0; ev<gNevents; ev++){
 347     if(skipEvt[ev]) continue;
 348     printf(" --- Processing event %d ---\n",ev);
 349
 350     gAlice->GetEvent(ev);
 351
 352     trk->cd();
 353
 354     // Tree with ITS tracks
 355     char tname[100];
 356     sprintf(tname,"TreeT_ITS_%d",ev);
 357
 358     TTree *tracktree=(TTree*)trk->Get(tname);
 359     if(!tracktree) continue;
 360     AliITStrackV2 *itstrack=new AliITStrackV2;
 361     tracktree->SetBranchAddress("tracks",&itstrack);
 362     Int_t nentr=(Int_t)tracktree->GetEntries();
 363
 364     // Tree for true track parameters
 365     char ttname[100];
 366     sprintf(ttname,"Tree_Ref_%d",ev);
 367     TTree *reftree = new TTree(ttname,"Tree with true track params");
 368     reftree->Branch("rectracks",&rectrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
 369
 370     for(Int_t i=0; i<nentr; i++) {
 371       tracktree->GetEvent(i);
 372       label = TMath::Abs(itstrack->GetLabel());
 373
 374       Part = (TParticle*)gAlice->Particle(label);
 375       rectrk.lab=label;
 376       rectrk.pdg=Part->GetPdgCode();
 377       rectrk.mumlab = Part->GetFirstMother();
 378       if(Part->GetFirstMother()>=0) {
 379         Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother());
 380         rectrk.mumpdg=Mum->GetPdgCode();
 381       } else {
 382         rectrk.mumpdg=-1;
 383       }
 384       rectrk.Vx=Part->Vx();
 385       rectrk.Vy=Part->Vy();
 386       rectrk.Vz=Part->Vz();
 387       rectrk.Px=Part->Px();
 388       rectrk.Py=Part->Py();
 389       rectrk.Pz=Part->Pz();
 390
 391       reftree->Fill();
 392     } // loop on tracks
 393
 394     out->cd();
 395     reftree->Write();
 396
 397     delete itstrack;
 398     delete reftree;
 399   } // loop on events
 400
 401   trk->Close();
 402   kin->Close();
 403   out->Close();
 404
 405   gBenchmark->Stop(name);
 406   gBenchmark->Show(name);
 407
 408
 409   return;
 410 }
 411 //-----------------------------------------------------------------------------
 412 void MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
 413
 414   printf("\n------------------------------------\n");
 415   printf("\nChecking for events to skip...\n");
 416
 417   Int_t evt,ncol;
 418
 419   FILE *f = fopen(evtsName,"r");
 420   while(1) {
 421     ncol = fscanf(f,"%d",&evt);
 422     if(ncol<1) break;
 423     skipEvt[evt] = 1;
 424     printf(" event %d will be skipped\n",evt);
 425   }
 426   fclose(f);
 427
 428   return;
 429 }
 430 //-----------------------------------------------------------------------------
 431 void PrimaryVertex(const Char_t *outName,Char_t vtxMode) {
 432
 433   printf("\n------------------------------------\n");
 434
 435   const Char_t *name="PrimaryVertex";
 436   printf("\n %s\n",name);
 437   gBenchmark->Start(name);
 438
 439   switch(vtxMode) {
 440   case 'H':
 441     printf(" ... from event header\n");
 442     VtxFromHeader(outName,kFALSE);
 443     break;
 444   case 'S':
 445     printf(" ... from event header + smearing\n");
 446     VtxFromHeader(outName,kTRUE);
 447     break;
 448   case 'P':
 449     printf(" ... z from pixels for pp/pA\n");
 450     ZvtxFromSPD(outName);
 451     break;
 452   case 'T':
 453     printf(" ... from tracks for pp/pA\n");
 454     VtxFromTracks(outName);
 455     break;
 456   case 'C':
 457     printf(" ... copied from Vtx4Tracking.root to AliITStracksV2.root\n");
 458     CopyVtx("Vtx4Tracking.root",outName);
 459     break;
 460   }
 461
 462   gBenchmark->Stop(name);
 463   gBenchmark->Show(name);
 464
 465   return;
 466 }
 467 //-----------------------------------------------------------------------------
 468 void TPCParamTracks(Int_t coll,Double_t Bfield) {
 469
 470   printf("\n------------------------------------\n");
 471
 472   const Char_t *name="TPCParamTracks";
 473   printf("\n %s\n",name);
 474   gBenchmark->Start(name);
 475
 476   TFile *outFile=TFile::Open("AliTPCtracksParam.root","recreate");
 477   TFile *inFile =TFile::Open("galice.root");
 478
 479   AliTPCtrackerParam tracker(coll,Bfield,gNevents);
 480   tracker.BuildTPCtracks(inFile,outFile);
 481
 482   delete gAlice; gAlice=0;
 483
 484   inFile->Close();
 485   outFile->Close();
 486
 487   gBenchmark->Stop(name);
 488   gBenchmark->Show(name);
 489
 490   return;
 491 }
 492 //-----------------------------------------------------------------------------
 493 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName) {
 494
 495     if(!gSystem->AccessPathName(logName,kFileExists)) {
 496       FILE *ifile = fopen(logName,"r");
 497       Int_t lEvt=0,nCol=1;
 498       while(nCol>0) {
 499         nCol = fscanf(ifile,"%d",&lEvt);
 500       }
 501       fclose(ifile);
 502       if(lEvt==gNevents) {
 503         printf(" All events already reconstructed\n");
 504         return 0;
 505       } else {
 506         FILE *ofile = fopen("evtsToSkip.dat","a");
 507         fprintf(ofile,"%d\n",lEvt);
 508         fclose(ofile);
 509       }
 510     } else {
 511       printf("File itstracking.log not found\n");
 512     }
 513
 514     return 1;
 515 }
 516 //-----------------------------------------------------------------------------
 517 void VtxFromHeader(const Char_t *outName,Bool_t smear) {
 518
 519   TDatime t;
 520   UInt_t seed = t.Get();
 521   gRandom->SetSeed(seed);
 522
 523   TFile *galice  = new TFile("galice.root");
 524   TFile *outFile = new TFile(outName,"update");
 525
 526   TDirectory *curdir;
 527   Double_t pos[3],sigma[3];
 528   if(smear) {
 529     sigma[0]=15.e-4;
 530     sigma[1]=15.e-4;
 531     sigma[2]=10.e-4;
 532   } else {
 533     sigma[0]=0.;
 534     sigma[1]=0.;
 535     sigma[2]=0.;
 536   }
 537   Char_t vname[20];
 538
 539   galice->cd();
 540
 541   for(Int_t ev=0; ev<gNevents; ev++){
 542     printf(" event %d\n",ev);
 543     sprintf(vname,"Vertex_%d",ev);
 544     TArrayF o = 0;
 545     o.Set(3);
 546     AliHeader* header = 0;
 547     TTree* treeE = (TTree*)gDirectory->Get("TE");
 548     treeE->SetBranchAddress("Header",&header);
 549     treeE->GetEntry(ev);
 550     AliGenEventHeader* genHeader = header->GenEventHeader();
 551     if(genHeader) {
 552       // get primary vertex position
 553       genHeader->PrimaryVertex(o);
 554       pos[0] = (Double_t)o[0];
 555       pos[1] = (Double_t)o[1];
 556       pos[2] = (Double_t)o[2];
 557       if(smear) {
 558         pos[0] = gRandom->Gaus(pos[0],sigma[0]);
 559         pos[1] = gRandom->Gaus(pos[1],sigma[1]);
 560         pos[2] = gRandom->Gaus(pos[2],sigma[2]);
 561       }
 562       // create AliITSVertex
 563       AliITSVertex *vertex = new AliITSVertex(pos,sigma,vname);
 564     } else {
 565       printf(" ! event header not found : setting vertex to (0,0,0) !");
 566       pos[0] = 0.;
 567       pos[1] = 0.;
 568       pos[2] = 0.;
 569       // create AliITSVertex
 570       AliITSVertex *vertex = new AliITSVertex(pos,sigma,vname);
 571     }
 572     delete header;
 573     // write AliITSVertex to file
 574     curdir = gDirectory;
 575     outFile->cd();
 576     if(smear) {
 577       vertex->SetTitle("vertex from header, smeared");
 578     } else {
 579       vertex->SetTitle("vertex from header");
 580     }
 581     vertex->Write();
 582     curdir->cd();
 583     vertex = 0;
 584   }
 585
 586   outFile->Close();
 587   galice->Close();
 588
 589   delete outFile;
 590   delete galice;
 591
 592   return;
 593 }
 594 //-----------------------------------------------------------------------------
 595 void VtxFromTracks(const Char_t *outName) {
 596
 597   // Open input and output files
 598   TFile *inFile  = new TFile("AliITStracksV2.root");
 599   TFile *outFile = new TFile(outName,"update");
 600
 601   // set AliRun object to 0
 602   if(gAlice) gAlice = 0;
 603
 604   // Create vertexer
 605   AliITSVertexerTracks *vertexer =
 606     new AliITSVertexerTracks(inFile,outFile,gBfieldValue);
 607   vertexer->SetFirstEvent(0);
 608   vertexer->SetLastEvent(gNevents-1);
 609   vertexer->SetDebug(0);
 610   vertexer->PrintStatus();
 611   // Find vertices
 612   vertexer->FindVertices();
 613
 614   delete vertexer;
 615
 616   inFile->Close();
 617   outFile->Close();
 618   delete inFile;
 619   delete outFile;
 620
 621   return;
 622 }
 623 //-----------------------------------------------------------------------------
 624 void ZvtxFromSPD(const Char_t *outName) {
 625
 626   // create fast RecPoints, which are used for vertex finding
 627   cerr<<"Fast AliITSRecPoint(s) !\n";
 628   gROOT->LoadMacro("$(ALICE_ROOT)/ITS/AliITSHits2FastRecPoints.C");
 629   AliITSHits2FastRecPoints(0,gNevents-1);
 630
 631   // delphi ---> azimuthal range to accept tracklets
 632   // window ---> window in Z around the peak of tracklets proj. in mm
 633   Float_t delphi=0.05;
 634   Float_t window=3.;
 635   Float_t initx=0.;
 636   Float_t inity=0.;
 637
 638   TFile *infile = new TFile("galice.root");
 639   TFile *outfile = new TFile(outName,"update");
 640
 641   AliITSVertexerPPZ *vertexer = new AliITSVertexerPPZ(infile,outfile,initx,inity);
 642   vertexer->SetFirstEvent(0);
 643   vertexer->SetLastEvent(gNevents-1);
 644   vertexer->SetDebug(0);
 645   vertexer->SetDiffPhiMax(delphi);
 646   vertexer->SetWindow(window);
 647   vertexer->PrintStatus();
 648   vertexer->FindVertices();
 649   delete vertexer;
 650   vertexer=0;
 651
 652   outfile->Close();
 653   infile->Close();
 654   delete infile;
 655   delete outfile;
 656
 657
 658   return;
 659 }
 660 //-----------------------------------------------------------------------------
 661
 662
 663
 664
 665
 666
 667
 668
 669