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) Fast points in ITS                                        *
   8  *             3) ITS cluster finding V2                                    *
   9  *             4) Determine z position of primary vertex                    *
  10  *                - read from event header in PbPb events                   *
  11  *                - determined using points on SPD in pp events (M.Masera)  *
  12  *             5) ITS track finding V2                                      *
  13  *                * in pp, redetermine the z position of the primary vertex *
  14  *                  using the reconstructed tracks                          *
  15  *             6) Create a reference file with simulation info (p,PDG...)   *
  16  *                                                                          *
  17  * If TString mode="A" all 6 steps are executed                             *
  18  * If TString mode="B" only steps 3-6 are executed                          *
  19  *                                                                          *
  20  * (Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it                    *
  21  *  from AliTPCtest.C & AliITStestV2.C by I.Belikov                         *
  22  ****************************************************************************/
  23 #if !defined(__CINT__) || defined(__MAKECINT__)
  24 //-- --- standard headers-------------
  25 #include <iostream.h>
  26 //--------Root headers ---------------
  27 #include <TFile.h>
  28 #include <TStopwatch.h>
  29 #include <TObject.h>
  30 #include <TSystem.h>
  31 #include <TH1.h>
  32 #include <TH2.h>
  33 #include <TProfile.h>
  34 #include <TTree.h>
  35 #include <TBranch.h>
  36 #include <TParticle.h>
  37 #include <TRandom.h>
  38 //----- AliRoot headers ---------------
  39 #include "alles.h"
  40 #include "AliRun.h"
  41 #include "AliHeader.h"
  42 #include "AliGenEventHeader.h"
  43 #include "AliMagF.h"
  44 #include "AliModule.h"
  45 #include "AliArrayI.h"
  46 #include "AliDigits.h"
  47 #include "AliITS.h"
  48 #include "AliTPC.h"
  49 #include "AliITSgeom.h"
  50 #include "AliITSRecPoint.h"
  51 #include "AliITSclusterV2.h"
  52 #include "AliITSsimulationFastPoints.h"
  53 #include "AliITStrackerV2.h"
  54 #include "AliKalmanTrack.h"
  55 #include "AliTPCtrackerParam.h"
  56 //-------------------------------------
  57 #endif
  58
  59 // structure for track references
  60 typedef struct {
  61   Int_t lab;
  62   Int_t pdg;
  63   Int_t mumlab;
  64   Int_t mumpdg;
  65   Float_t Vx,Vy,Vz;
  66   Float_t Px,Py,Pz;
  67 } RECTRACK;
  68
  69 //===== Functions definition =================================================
  70
  71 Int_t TPCParamTracks(const Char_t *galiceName,const Char_t *outName,
  72                      const Int_t coll,const Double_t Bfield,Int_t n);
  73
  74 Int_t ITSHits2FastRecPoints(const Char_t *galiceName,Int_t n);
  75
  76 Int_t ITSFindClustersV2(const Char_t *galiceName,const Char_t *outName,
  77                         Int_t *skipEvt,Int_t n);
  78
  79 Int_t ZvtxFromHeader(const Char_t *galiceName,Int_t n);
  80
  81 Int_t ZvtxFromSPD(const Char_t *galiceName,Int_t n);
  82
  83 Int_t ZvtxFromTracks(const Char_t *trkame,Int_t n);
  84
  85 Int_t ZvtxFastpp(const Char_t *galiceName,Int_t n);
  86
  87 void EvalZ(TH1F *hist,Int_t sepa,Float_t &av,Float_t &sig,Int_t ncoinc,
  88            TArrayF *zval,ofstream *deb);
  89
  90 Bool_t VtxTrack(Double_t pt,Double_t d0rphi);
  91
  92 Double_t d0zRes(Double_t pt);
  93
  94 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName,Int_t n);
  95
  96 Int_t MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt);
  97
  98 Int_t ITSFindTracksV2(const Char_t *galiceName,const Char_t *inName,
  99                       const Char_t *inName2,const Char_t *outName,
 100                       Int_t *skipEvt,TString vtxMode,Int_t n);
 101
 102 Int_t ITSMakeRefFile(const Char_t *galiceName,const Char_t *inName,
 103                      const Char_t *outName,Int_t *skipEvt,Int_t n);
 104
 105 void WriteVtx(const Char_t *name,Double_t *zvtx,Int_t n);
 106
 107 void ReadVtx(const Char_t *name,Double_t *zvtx,Int_t n);
 108
 109 //=============================================================================
 110
 111 void AliBarrelRec_TPCparam(TString mode="A",Int_t n=5000) {
 112
 113   const Char_t *name=" AliBarrelRec_TPCparam";
 114   cerr<<'\n'<<name<<"...\n";
 115   gBenchmark->Start(name);
 116
 117
 118   // filenames
 119   const Char_t *galiceName  = "galice.root";
 120   const Char_t *TPCtrkNameS = "AliTPCtracksParam.root";
 121   const Char_t *ITSclsName  = "AliITSclustersV2.root";
 122   const Char_t *ITStrkName  = "AliITStracksV2.root";
 123   const Char_t *ITStrkName2 = "AliITStracksV2_2.root";
 124   const Char_t *ITSrefName  = "ITStracksRefFile.root";
 125
 126   // set here the code for the type of collision (needed for TPC tracking
 127   // parameterization). available collisions:
 128   //
 129   // coll = 0 ->   PbPb6000 (HIJING with b<2fm)
 130   // coll = 1 ->   pp
 131   const Int_t    collcode    = 1;
 132   const Bool_t   slowVtx     = kFALSE;
 133   const Bool_t   retrack     = kFALSE;
 134   // set here the value of the magnetic field
 135   const Double_t BfieldValue = 0.4;
 136
 137   // set conversion constant for kalman tracks
 138   AliKalmanTrack::SetConvConst(100/0.299792458/BfieldValue);
 139
 140   Bool_t clust = kTRUE;
 141
 142   //################ SWITCH #############################################
 143   switch (*(mode.Data())) {
 144     //############## case A #############################################
 145   case 'A':
 146     // ********** Build TPC tracks with parameterization *********** //
 147     TPCParamTracks(galiceName,TPCtrkNameS,collcode,BfieldValue,n);
 148
 149     // ********** ITS RecPoints ************************************ //
 150     ITSHits2FastRecPoints(galiceName,n);
 151
 152     break;
 153     //############## case B #############################################
 154   case 'B':
 155     cerr<<"       ---> only tracking in ITS <---"<<endl;
 156
 157     if(!UpdateEvtsToSkip("itstracking.log","evtsToSkip.dat",n)) return;
 158
 159     if(!UpdateEvtsToSkip("itsclusters.log","evtsToSkip.dat",n)) clust=kFALSE;
 160
 161     break;
 162     //############## END SWITCH #########################################
 163   }
 164
 165   // Mark events that have to be skipped (if any)
 166   Int_t *skipEvt = new Int_t[n];
 167   for(Int_t i=0;i<n;i++) skipEvt[i] = 0;
 168   if(!gSystem->AccessPathName("evtsToSkip.dat",kFileExists)) {
 169     MarkEvtsToSkip("evtsToSkip.dat",skipEvt);
 170   }
 171
 172   // ********** Find ITS clusters ******************************** //
 173   if(clust) ITSFindClustersV2(galiceName,ITSclsName,skipEvt,n);
 174
 175   // ********** Tracking in ITS ********************************** //
 176   TString vtxMode;
 177   switch (collcode) {
 178   case 0:   // ***** Pb-Pb *****
 179     ZvtxFromHeader(galiceName,n);
 180     vtxMode="H";
 181     ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,
 182                     skipEvt,vtxMode,n);
 183     break;
 184   case 1:   // *****  pp   *****
 185     if(slowVtx) {
 186       ZvtxFromSPD(galiceName,n);
 187       vtxMode="P";
 188       ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,
 189                       skipEvt,vtxMode,n);
 190       ZvtxFromTracks(ITStrkName,n);
 191       if(retrack) {
 192         vtxMode="T";
 193         ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName2,
 194                         skipEvt,vtxMode,n);
 195       }
 196     } else {
 197       ZvtxFastpp(galiceName,n);
 198       vtxMode="F";
 199       ITSFindTracksV2(galiceName,TPCtrkNameS,ITSclsName,ITStrkName,
 200                       skipEvt,vtxMode,n);
 201     }
 202     break;
 203   }
 204
 205
 206   // ********** Make ITS tracks reference file ******************* //
 207   ITSMakeRefFile(galiceName,ITStrkName,ITSrefName,skipEvt,n);
 208
 209   delete [] skipEvt;
 210
 211
 212
 213   gBenchmark->Stop(name);
 214   gBenchmark->Show(name);
 215
 216   return;
 217 }
 218 //-----------------------------------------------------------------------------
 219 Int_t UpdateEvtsToSkip(const Char_t *logName,const Char_t *evtsName,Int_t n) {
 220
 221     if(!gSystem->AccessPathName(logName,kFileExists)) {
 222       FILE *ifile = fopen(logName,"r");
 223       Int_t lEvt=0,nCol=1;
 224       while(nCol>0) {
 225         nCol = fscanf(ifile,"%d",&lEvt);
 226       }
 227       fclose(ifile);
 228       if(lEvt==n) {
 229         cerr<<" All events already reconstructed "<<endl;
 230         return 0;
 231       } else {
 232         FILE *ofile = fopen("evtsToSkip.dat","a");
 233         fprintf(ofile,"%d\n",lEvt);
 234         fclose(ofile);
 235       }
 236     } else {
 237       cerr<<"File itstracking.log not found"<<endl;
 238     }
 239
 240     return 1;
 241 }
 242 //-----------------------------------------------------------------------------
 243 Int_t MarkEvtsToSkip(const Char_t *evtsName,Int_t *skipEvt) {
 244
 245   cerr<<"\n*******************************************************************\n";
 246   cerr<<"\nChecking for events to skip...\n";
 247
 248   Int_t evt,ncol;
 249
 250   FILE *f = fopen(evtsName,"r");
 251   while(1) {
 252     ncol = fscanf(f,"%d",&evt);
 253     if(ncol<1) break;
 254     skipEvt[evt] = 1;
 255     cerr<<" ev. "<<evt<<" will be skipped\n";
 256   }
 257   fclose(f);
 258
 259   return 0;
 260 }
 261 //-----------------------------------------------------------------------------
 262 Int_t TPCParamTracks(const Char_t *galiceName,const Char_t *outName,
 263                      const Int_t coll,const Double_t Bfield,Int_t n) {
 264
 265   cerr<<"\n*******************************************************************\n";
 266
 267   const Char_t *name="TPCParamTracks";
 268   cerr<<'\n'<<name<<"...\n";
 269   gBenchmark->Start(name);
 270
 271   TFile *outFile=TFile::Open(outName,"recreate");
 272   TFile *inFile =TFile::Open(galiceName);
 273
 274   AliTPCtrackerParam tracker(coll,Bfield,n);
 275   tracker.BuildTPCtracks(inFile,outFile);
 276
 277   delete gAlice; gAlice=0;
 278
 279   inFile->Close();
 280   outFile->Close();
 281
 282   gBenchmark->Stop(name);
 283   gBenchmark->Show(name);
 284
 285   return 0;
 286 }
 287 //-----------------------------------------------------------------------------
 288 Int_t ITSHits2FastRecPoints(const Char_t *galiceName,Int_t n) {
 289
 290   cerr<<"\n*******************************************************************\n";
 291
 292   const Char_t *name="ITSHits2FastRecPoints";
 293   cerr<<'\n'<<name<<"...\n";
 294   gBenchmark->Start(name);
 295
 296   Int_t nsignal=25;
 297   Int_t size=-1;
 298
 299   // Connect the Root Galice file containing Geometry, Kine and Hits
 300   TFile *file = TFile::Open(galiceName,"UPDATE");
 301
 302   gAlice = (AliRun*)file->Get("gAlice");
 303
 304
 305   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 306   if(!ITS) return 1;
 307
 308   // Set the simulation model
 309   for(Int_t i=0;i<3;i++) {
 310     ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
 311   }
 312
 313
 314   //
 315   // Event Loop
 316   //
 317   for(Int_t ev=0; ev<n; ev++) {
 318     cerr<<" --- Processing event "<<ev<<" ---"<<endl;
 319     Int_t nparticles = gAlice->GetEvent(ev);
 320     cerr<<"Number of particles: "<<nparticles<<endl;
 321     gAlice->SetEvent(ev);
 322     if(!gAlice->TreeR()) gAlice->MakeTree("R");
 323     if(!gAlice->TreeR()->GetBranch("ITSRecPointsF")) {
 324       ITS->MakeBranch("RF");
 325       if(nparticles <= 0) return 1;
 326
 327       Int_t bgr_ev=Int_t(ev/nsignal);
 328       //printf("bgr_ev %d\n",bgr_ev);
 329       ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
 330     }
 331   } // event loop
 332
 333   delete gAlice; gAlice=0;
 334   file->Close();
 335
 336   gBenchmark->Stop(name);
 337   gBenchmark->Show(name);
 338
 339   return 0;
 340 }
 341 //-----------------------------------------------------------------------------
 342 Int_t ITSFindClustersV2(const Char_t *galiceName,const Char_t *outName,
 343                         Int_t *skipEvt,Int_t n) {
 344   //
 345   // This function converts AliITSRecPoint(s) to AliITSclusterV2
 346   //
 347   cerr<<"\n*******************************************************************\n";
 348
 349   const Char_t *name="ITSFindClustersV2";
 350   cerr<<'\n'<<name<<"...\n";
 351   gBenchmark->Start(name);
 352
 353
 354   TFile *inFile=TFile::Open(galiceName);
 355   if(!inFile->IsOpen()) { cerr<<"Can't open galice.root !\n"; return 1; }
 356
 357   if(!(gAlice=(AliRun*)inFile->Get("gAlice"))) {
 358     cerr<<"Can't find gAlice !\n";
 359     return 1;
 360   }
 361   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 362   if(!ITS) { cerr<<"Can't find the ITS !\n"; return 1; }
 363   AliITSgeom *geom=ITS->GetITSgeom();
 364
 365   TFile *outFile = TFile::Open(outName,"recreate");
 366   geom->Write();
 367
 368   // open logfile for done events
 369   FILE *logfile = fopen("itsclusters.log","w");
 370
 371   // loop on events
 372   for(Int_t ev=0; ev<n; ev++) {
 373     // write to logfile of begun events
 374     fprintf(logfile,"%d\n",ev);
 375
 376     if(skipEvt[ev]) continue;
 377
 378     cerr<<"--- Processing event "<<ev<<" ---"<<endl;
 379     inFile->cd();
 380     gAlice->GetEvent(ev);
 381
 382     TClonesArray *clusters = new TClonesArray("AliITSclusterV2",10000);
 383     Char_t cname[100];
 384     sprintf(cname,"TreeC_ITS_%d",ev);
 385     TTree *cTree = new TTree(cname,"ITS clusters");
 386     cTree->Branch("Clusters",&clusters);
 387
 388     TTree *pTree=gAlice->TreeR();
 389     if(!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
 390     TBranch *branch = pTree->GetBranch("ITSRecPointsF");
 391     if(!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1; }
 392     TClonesArray *points = new TClonesArray("AliITSRecPoint",10000);
 393     branch->SetAddress(&points);
 394
 395     TClonesArray &cl=*clusters;
 396     Int_t nclusters=0;
 397     Int_t nentr=(Int_t)branch->GetEntries();
 398
 399     cerr<<"Number of entries in TreeR_RF: "<<nentr<<endl;
 400
 401     //    Float_t *lp  = new Float_t[5];
 402     //    Int_t   *lab = new Int_t[6];
 403     Float_t lp[5];
 404     Int_t   lab[6];
 405
 406     for(Int_t i=0; i<nentr; i++) {
 407       points->Clear();
 408       branch->GetEvent(i);
 409       Int_t ncl=points->GetEntriesFast(); if(ncl==0){cTree->Fill();continue;}
 410       Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
 411       Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
 412       Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
 413       Double_t yshift = x*rot[0] + y*rot[1];
 414       Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
 415       nclusters+=ncl;
 416
 417       Float_t kmip=1; // ADC->mip normalization factor for the SDD and SSD
 418       if(lay==4 || lay==3) kmip=280.;
 419       if(lay==6 || lay==5) kmip= 38.;
 420
 421       for(Int_t j=0; j<ncl; j++) {
 422         AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
 423         lp[0]=-p->GetX()-yshift; if(lay==1) lp[0]=-lp[0];
 424         lp[1]=p->GetZ()+zshift;
 425         lp[2]=p->GetSigmaX2();
 426         lp[3]=p->GetSigmaZ2();
 427         lp[4]=p->GetQ(); lp[4]/=kmip;
 428         lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
 429         lab[3]=ndet;
 430
 431         Int_t label=lab[0];
 432         if(label>=0) {
 433           TParticle *part=(TParticle*)gAlice->Particle(label);
 434           label=-3;
 435           while (part->P() < 0.005) {
 436             Int_t m=part->GetFirstMother();
 437             if(m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
 438             label=m;
 439             part=(TParticle*)gAlice->Particle(label);
 440           }
 441           if      (lab[1]<0) lab[1]=label;
 442           else if(lab[2]<0) lab[2]=label;
 443           else cerr<<"No empty labels !\n";
 444         }
 445         new(cl[j]) AliITSclusterV2(lab,lp);
 446       }
 447       cTree->Fill(); clusters->Delete();
 448       points->Delete();
 449     }
 450
 451     outFile->cd();
 452     cTree->Write();
 453
 454     cerr<<"Number of clusters: "<<nclusters<<endl;
 455     //delete [] lp;
 456     //delete [] lab;
 457     delete cTree; delete clusters; delete points;
 458
 459   } // end loop on events
 460
 461   fprintf(logfile,"%d\n",n); //this means all evts are successfully completed
 462   fclose(logfile);
 463
 464   delete gAlice; gAlice=0;
 465
 466   inFile->Close();
 467   outFile->Close();
 468
 469   gBenchmark->Stop(name);
 470   gBenchmark->Show(name);
 471
 472    return 0;
 473 }
 474 //-----------------------------------------------------------------------------
 475 Int_t ZvtxFromHeader(const Char_t *galiceName,Int_t n) {
 476
 477   cerr<<"\n*******************************************************************\n";
 478
 479   const Char_t *name="ZvtxFromHeader";
 480   cerr<<'\n'<<name<<"...\n";
 481   gBenchmark->Start(name);
 482
 483   TFile *galice = TFile::Open(galiceName);
 484
 485   Double_t *zvtx    = new Double_t[n];
 486
 487   for(Int_t ev=0; ev<n; ev++){
 488     cerr<<" --- Processing event "<<ev<<" ---"<<endl;
 489
 490     TArrayF o = 0;
 491     o.Set(3);
 492     AliHeader* header = 0;
 493     TTree* treeE = (TTree*)gDirectory->Get("TE");
 494     treeE->SetBranchAddress("Header",&header);
 495     treeE->GetEntry(ev);
 496     AliGenEventHeader* genHeader = header->GenEventHeader();
 497     if(genHeader) {
 498       // get primary vertex position
 499       genHeader->PrimaryVertex(o);
 500       // set position of primary vertex
 501       zvtx[ev] = (Double_t)o[2];
 502     } else {
 503       cerr<<" ! event header not found : setting z vertex to 0 !"<<endl;
 504       zvtx[ev] = 0.;
 505     }
 506     delete header;
 507   }
 508
 509   galice->Close();
 510
 511   // Write vertices to file
 512   WriteVtx("vtxHeader.root",zvtx,n);
 513
 514   delete [] zvtx;
 515
 516   gBenchmark->Stop(name);
 517   gBenchmark->Show(name);
 518
 519   return 0;
 520 }
 521 //-----------------------------------------------------------------------------
 522 Int_t ZvtxFastpp(const Char_t *galiceName,Int_t n) {
 523
 524   cerr<<"\n*******************************************************************\n";
 525
 526   const Char_t *name="ZvtxFastpp";
 527   cerr<<'\n'<<name<<"...\n";
 528   gBenchmark->Start(name);
 529
 530   TFile *galice = TFile::Open(galiceName);
 531   AliRun *gAlice = (AliRun*)galice->Get("gAlice");
 532
 533   Double_t *zvtx    = new Double_t[n];
 534
 535   TParticle *part;
 536   Int_t      nPart,b;
 537   Double_t dNchdy,E,theta,eta,zvtxTrue,sigmaz,probVtx;
 538
 539   for(Int_t ev=0; ev<n; ev++) {
 540     cerr<<" --- Processing event "<<ev<<" ---"<<endl;
 541
 542     TArrayF o = 0;
 543     o.Set(3);
 544     AliHeader* header = 0;
 545     TTree* treeE = (TTree*)galice->Get("TE");
 546     treeE->SetBranchAddress("Header",&header);
 547     treeE->GetEntry(ev);
 548     AliGenEventHeader* genHeader = header->GenEventHeader();
 549     if(genHeader) {
 550       // get primary vertex position
 551       genHeader->PrimaryVertex(o);
 552       // set position of primary vertex
 553       zvtxTrue = (Double_t)o[2];
 554     } else {
 555       cerr<<" ! event header not found : setting z vertex to 0 !"<<endl;
 556       zvtxTrue = 0.;
 557     }
 558     delete header;
 559
 560
 561     // calculate dNch/dy pf the event (charged primaries in |eta|<0.5)
 562     nPart = (Int_t)gAlice->GetEvent(ev);
 563     dNchdy = 0.;
 564     for(b=0; b<nPart; b++) {
 565       part = (TParticle*)gAlice->Particle(b);
 566       if(TMath::Abs(part->GetPdgCode())<10) continue; // reject partons
 567       if(TMath::Abs(part->Vx()*part->Vx()+part->Vy()*part->Vy())>0.01) continue; // reject secondaries
 568       E  = part->Energy();
 569       if(E>6900.) continue; // reject incoming protons
 570       theta = part->Theta();
 571       if(theta<.1 || theta>3.) continue;
 572       eta = -TMath::Log(TMath::Tan(theta/2.));
 573       if(TMath::Abs(eta)>0.5) continue; // count particles in |eta|<0.5
 574       dNchdy+=1.;
 575     }
 576
 577     // get sigma(z) corresponding to the mult of this event
 578     if(dNchdy>1.)  {
 579       sigmaz   = 0.0417/TMath::Sqrt(dNchdy);
 580     } else {
 581       sigmaz = 0.0500;
 582     }
 583     // smear the original position of the primary vertex
 584     zvtx[ev] = gRandom->Gaus(zvtxTrue,sigmaz);
 585
 586     // compute the probability that the vertex is found
 587     probVtx = 1.;
 588     if(dNchdy<24.) probVtx = 0.85+0.00673*dNchdy;
 589     if(dNchdy<14.) probVtx = 0.85;
 590
 591     if(gRandom->Rndm()>probVtx) zvtx[ev] = -1000.;// no zvtx for this event
 592
 593   }
 594   galice->Close();
 595
 596   // Write vertices to file
 597   WriteVtx("vtxFastpp.root",zvtx,n);
 598
 599   delete [] zvtx;
 600   //delete gAlice; gAlice=0;
 601
 602
 603   gBenchmark->Stop(name);
 604   gBenchmark->Show(name);
 605
 606   return 0;
 607 }
 608 //-----------------------------------------------------------------------------
 609 Int_t ZvtxFromSPD(const Char_t *galiceName,Int_t n) {
 610
 611
 612   Double_t *zvtx    = new Double_t[n];
 613
 614   Bool_t debug = kFALSE;
 615   ofstream fildeb;
 616   if(debug)fildeb.open("FindZV.out");
 617   ofstream filou;
 618   //filou.open("zcoor.dat");
 619   const Float_t kPi=TMath::Pi();
 620   const Int_t kFirstL1=0;
 621   const Int_t kLastL1=79;
 622   const Int_t kFirstL2=80;
 623   const Int_t kLastL2=239;
 624   //  const Float_t kDiffPhiMax=0.0005;
 625   const Float_t kDiffPhiMax=0.01;
 626   const Float_t kphl=0.;
 627   const Float_t kphh=kPi/18.;
 628   TDirectory *currdir = gDirectory;
 629   TH2F *hz2z1 = 0;
 630   TH1F *zvdis = 0;
 631   TH2F *dpvsz = 0;
 632   TProfile *scoc = new TProfile("scoc","Zgen - Zmeas vs occ lay 1",5,10.,60.);
 633   TProfile *prof = new TProfile("prof","z2 vs z1",50,-15.,15.);
 634   TH1F *phdiff = new TH1F("phdiff","#Delta #phi",100,kphl,kphh);
 635   TH1F *phdifsame = new TH1F("phdifsame","#Delta #phi same track",100,kphl,kphh);
 636   if(gAlice){delete gAlice; gAlice = 0;}
 637   TFile *file = TFile::Open(galiceName);
 638   gAlice = (AliRun*)file->Get("gAlice");
 639
 640   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 641   AliITSgeom *geom = ITS->GetITSgeom();
 642   TTree *TR=0;
 643   Float_t lc[3]; for(Int_t ii=0; ii<3; ii++) lc[ii]=0.;
 644   Float_t gc[3]; for(Int_t ii=0; ii<3; ii++) gc[ii]=0.;
 645   Float_t lc2[3]; for(Int_t ii=0; ii<3; ii++) lc2[ii]=0.;
 646   Float_t gc2[3]; for(Int_t ii=0; ii<3; ii++) gc2[ii]=0.;
 647   char name[30];
 648   Float_t avesca=0;
 649   Float_t avesca2=0;
 650   Int_t N=0;
 651   for(Int_t event=0; event<n; event++){
 652     sprintf(name,"event_%d",event);
 653     hz2z1=new TH2F(name,"z2 vs z1",50,-15.,15.,50,-15.,15.);
 654     hz2z1->SetDirectory(currdir);
 655     gAlice->GetEvent(event);
 656     TParticle * part = gAlice->Particle(0);
 657     Float_t oriz=part->Vz();
 658     cout<<"\n ==============================================================\n";
 659     cout<<"   Processing event "<<event<<" "<<oriz<<endl;
 660     cout<<" ==============================================================\n";
 661     if(debug){
 662     fildeb<<"\n ==============================================================\n";
 663     fildeb<<"   Processing event "<<event<<" "<<oriz<<endl;
 664     fildeb<<" ==============================================================\n";
 665     }
 666     file->cd();
 667     TR = gAlice->TreeR();
 668     if(!TR) cerr<<"TreeR not found\n";
 669     TBranch *branch=TR->GetBranch("ITSRecPointsF");
 670     if(!branch) cerr<<"Branch ITSRecPointsF not found\n";
 671     TClonesArray *points  = new TClonesArray("AliITSRecPoint",10000);
 672     branch->SetAddress(&points);
 673     Int_t nentr=(Int_t)branch->GetEntries();
 674     Float_t zave=0;
 675     Float_t rmszav=0;
 676     Float_t zave2=0;
 677     Int_t firipixe=0;
 678     cout<<endl;
 679     for(Int_t module= kFirstL1; module<=kLastL1;module++){
 680       points->Clear();
 681       branch->GetEvent(module);
 682       Int_t nrecp1 = points->GetEntriesFast();
 683       //cerr<<nrecp1<<endl;
 684       for(Int_t i=0; i<nrecp1;i++){
 685         AliITSRecPoint *current = (AliITSRecPoint*)points->UncheckedAt(i);
 686         lc[0]=current->GetX();
 687         lc[2]=current->GetZ();
 688         geom->LtoG(module,lc,gc);
 689         zave+=gc[2];
 690         zave2+=gc[2]*gc[2];
 691         firipixe++;
 692       }
 693       points->Delete();
 694     }
 695     if(firipixe>1){
 696       rmszav=TMath::Sqrt(zave2/(firipixe-1)-zave*zave/firipixe/(firipixe-1));
 697       zave=zave/firipixe;
 698       cout<<"Z aver first layer = "<<zave<<" RMS= "<<rmszav<<" pixels = "<<firipixe<<endl;
 699       if(debug)fildeb<<"Z aver first layer = "<<zave<<" RMS= "<<rmszav<<" pixels = "<<firipixe<<endl;
 700     }
 701     else {
 702       cout<<"No rec points on first layer for this event\n";
 703       if(debug)fildeb<<"No rec points on first layer for this event\n";
 704     }
 705     Float_t zlim1=zave-rmszav;
 706     Float_t zlim2=zave+rmszav;
 707     Int_t sepa=(Int_t)((zlim2-zlim1)*10.+1.);
 708     zlim2=zlim1 + sepa/10.;
 709     sprintf(name,"pz_ev_%d",event);
 710     dpvsz=new TH2F(name,"#Delta #phi vs Z",sepa,zlim1,zlim2,100,0.,0.03);
 711     dpvsz->SetDirectory(currdir);
 712     sprintf(name,"z_ev_%d",event);
 713     zvdis=new TH1F(name,"zv distr",sepa,zlim1,zlim2);
 714     zvdis->SetDirectory(currdir);
 715     cout<<"Z limits: "<<zlim1<<" "<<zlim2<<endl;
 716     if(debug)fildeb<<"Z limits: "<<zlim1<<" "<<zlim2<<endl;
 717     Int_t sizarr=100;
 718     TArrayF *zval = new TArrayF(sizarr);
 719     Int_t ncoinc=0;
 720     for(Int_t module= kFirstL1; module<=kLastL1;module++){
 721       //cout<<"processing module   "<<module<<"                  \r";
 722       branch->GetEvent(module);
 723       Int_t nrecp1 = points->GetEntriesFast();
 724       TObjArray *poiL1 = new TObjArray(nrecp1);
 725       for(Int_t i=0; i<nrecp1;i++)poiL1->AddAt(points->UncheckedAt(i),i);
 726       //ITS->ResetRecPoints();
 727       points->Delete();
 728       for(Int_t i=0; i<nrecp1;i++){
 729         AliITSRecPoint *current = (AliITSRecPoint*)poiL1->UncheckedAt(i);
 730         lc[0]=current->GetX();
 731         lc[2]=current->GetZ();
 732         geom->LtoG(module,lc,gc);
 733         Float_t r1=TMath::Sqrt(gc[0]*gc[0]+gc[1]*gc[1]);
 734         Float_t phi1 = TMath::ATan2(gc[1],gc[0]);
 735         if(phi1<0)phi1=2*kPi+phi1;
 736         Int_t lab1 = current->GetLabel(0);
 737         Float_t phi1d=phi1*180./kPi;
 738         //cerr<<"module "<<module<<" "<<gc[0]<<" "<<gc[1]<<" "<<gc[2]<<" "<<phi1d<<" "<<lab1<<"     \n";
 739         for(Int_t modul2=kFirstL2; modul2<=kLastL2; modul2++){
 740           branch->GetEvent(modul2);
 741           Int_t nrecp2 = points->GetEntriesFast();
 742           for(Int_t j=0; j<nrecp2;j++){
 743             AliITSRecPoint *recp = (AliITSRecPoint*)points->UncheckedAt(j);
 744             lc2[0]=recp->GetX();
 745             lc2[2]=recp->GetZ();
 746             geom->LtoG(modul2,lc2,gc2);
 747             Float_t r2=TMath::Sqrt(gc2[0]*gc2[0]+gc2[1]*gc2[1]);
 748             Float_t zr0=(r2*gc[2]-r1*gc2[2])/(r2-r1);
 749             Int_t lab2 = recp->GetLabel(0);
 750             Float_t phi2 = TMath::ATan2(gc2[1],gc2[0]);
 751             if(phi2<0)phi2=2.*kPi+phi2;
 752             Float_t diff = TMath::Abs(phi2-phi1);
 753             if(diff>kPi)diff=2.*kPi-diff;
 754             if(lab1==lab2)phdifsame->Fill(diff);
 755             if(zr0>zlim1 && zr0<zlim2){
 756               phdiff->Fill(diff);
 757               dpvsz->Fill(zr0,diff);
 758               if(diff<kDiffPhiMax ){
 759                 zvdis->Fill(zr0);
 760                 zval->AddAt(zr0,ncoinc);
 761                 ncoinc++;
 762                 if(ncoinc==(sizarr-1)){
 763                   sizarr+=100;
 764                   zval->Set(sizarr);
 765                 }
 766                 Float_t phi2d=phi2*180./kPi;
 767                 Float_t diffd=diff*180./kPi;
 768                 //cerr<<"module2 "<<modul2<<" "<<gc2[0]<<" "<<gc2[1]<<" "<<gc2[2]<<" "<<phi2d<<" "<<diffd<<" "<<lab2<<"     \n";
 769                 //                cout<<"zr0= "<<zr0<<endl;
 770                 hz2z1->Fill(gc[2],gc2[2]);
 771                 prof->Fill(gc[2]-oriz,gc2[2]-oriz);
 772               }
 773             }
 774           }
 775           points->Delete();
 776         }
 777       }
 778       delete poiL1;
 779     }         // loop on modules
 780     cout<<endl;
 781     Float_t zcmp=0;
 782     Float_t zsig=0;
 783     EvalZ(zvdis,sepa,zcmp,zsig,ncoinc,zval,&fildeb);
 784     if(zcmp!=-100 && zsig!=-100){
 785       N++;
 786       Float_t scarto = (oriz-zcmp)*10000.;
 787       Float_t ascarto=TMath::Abs(scarto);
 788       scoc->Fill(firipixe,ascarto);
 789       avesca+=ascarto;
 790       avesca2+=scarto*scarto;
 791       if(debug)fildeb<<"Mean value: "<<zcmp<<"  +/-"<<zsig<<"  Zgen- Zmeas (micron)= "<<scarto<<endl;
 792       cout<<"Mean value: "<<zcmp<<"  RMS: "<<zsig<<"  Zgen- Zmeas (micron)= "<<scarto<<endl;
 793       //filou<<event<<" "<<zcmp<<" "<<zsig<<" "<<oriz<<" "<<scarto<<endl;
 794       zvtx[event] = (Double_t)zcmp;
 795     }
 796     else {
 797       cout<<"Not enough points in event "<<event<<endl;
 798       if(debug)fildeb<<"Not enough points\n";
 799       //filou<<event<<" "<<-1000.<<" "<<0.<<" "<<oriz<<" "<<0.<<endl;
 800       zvtx[event] = -1000.;
 801     }
 802     delete zval;
 803   }           // loop on events
 804   file->Close();
 805   //hz2z1->Draw("box");
 806
 807   // Write vertices to file
 808   WriteVtx("vtxSPD.root",zvtx,n);
 809   delete [] zvtx;
 810
 811   if(N>1) {
 812     avesca2=TMath::Sqrt((avesca2/(N-1)-avesca*avesca/N/(N-1))/N);
 813     avesca=avesca/N;
 814   } else {
 815     avesca2 = 0;
 816   }
 817   cout<<"\n \n ==========================================================\n";
 818   cout<<"Number of events with vertex estimate  "<<N<<endl;
 819   cout<<"Average of the (abs) Zgen-Zmeas   :   "<<avesca;
 820   cout<<" +/- "<<avesca2<<" micron"<<endl;
 821   if(debug){
 822     fildeb<<"\n \n ==========================================================\n";
 823     fildeb<<"Number of events with vertex estimate  "<<N<<endl;
 824     fildeb<<"Average of the (abs) Zgen-Zmeas   :   "<<avesca;
 825     fildeb<<" +/- "<<avesca2<<endl;
 826   }
 827   return 0;
 828 }
 829
 830
 831 //-----------------------------------------------------------------------------
 832 void EvalZ(TH1F *hist,Int_t sepa,Float_t &av, Float_t &sig,Int_t ncoinc,
 833            TArrayF *zval,ofstream *deb) {
 834
 835   av=0;
 836   sig=0;
 837   Int_t N=0;
 838   Int_t NbinNotZero=0;
 839   for(Int_t i=1;i<=sepa;i++){
 840     Float_t cont=hist->GetBinContent(i);
 841     av+=cont;
 842     sig+=cont*cont;
 843     N++;
 844     if(cont!=0)NbinNotZero++;
 845   }
 846   if(NbinNotZero==0){av=-100; sig=-100; return;}
 847   if(NbinNotZero==1){sig=-100; return;}
 848   sig=TMath::Sqrt(sig/(N-1)-av*av/N/(N-1));
 849   av=av/N;
 850   Float_t val1=hist->GetBinLowEdge(sepa);
 851   Float_t val2=hist->GetBinLowEdge(1);
 852   for(Int_t i=1;i<=sepa;i++){
 853     Float_t cont=hist->GetBinContent(i);
 854     if(cont>(av+sig*2.)){
 855       Float_t curz=hist->GetBinLowEdge(i);
 856       if(curz<val1)val1=curz;
 857       if(curz>val2)val2=curz;
 858     }
 859   }
 860   val2+=hist->GetBinWidth(1);
 861   cout<<"Values for Z finding: "<<val1<<" "<<val2<<endl;
 862   if(deb)(*deb)<<"Values for Z finding: "<<val1<<" "<<val2<<endl;
 863   av=0;
 864   sig=0;
 865   N=0;
 866   for(Int_t i=0; i<ncoinc; i++){
 867     Float_t z=zval->At(i);
 868     if(z<val1)continue;
 869     if(z>val2)continue;
 870     av+=z;
 871     sig+=z*z;
 872     N++;
 873   }
 874   if(N<=1){av=-100; sig=-100; return;}
 875   sig=TMath::Sqrt((sig/(N-1)-av*av/N/(N-1))/N);
 876   av=av/N;
 877   return;
 878 }
 879 //-----------------------------------------------------------------------------
 880 Int_t ZvtxFromTracks(const Char_t *trkName,Int_t n) {
 881
 882   cerr<<"\n*******************************************************************\n";
 883
 884   const Char_t *name="ZvtxFromTracks";
 885   cerr<<'\n'<<name<<"...\n";
 886   gBenchmark->Start(name);
 887
 888   Double_t *zvtx = new Double_t[n];
 889   Double_t *zvtxSPD = new Double_t[n];
 890   ReadVtx("vtxSPD.root",zvtxSPD,n);
 891
 892   TFile *itstrk = TFile::Open(trkName);
 893
 894
 895   Int_t    nTrks,i,nUsedTrks;
 896   Double_t pt,d0rphi,d0z;
 897   Double_t zvtxTrks,ezvtxTrks,sumWeights;
 898
 899   for(Int_t ev=0; ev<n; ev++){
 900     cerr<<" --- Processing event "<<ev<<" ---"<<endl;
 901
 902     // Tree with ITS tracks
 903     Char_t tname[100];
 904     sprintf(tname,"TreeT_ITS_%d",ev);
 905
 906     TTree *tracktree=(TTree*)itstrk->Get(tname);
 907     if(!tracktree) continue;
 908     AliITStrackV2 *itstrack=new AliITStrackV2;
 909     tracktree->SetBranchAddress("tracks",&itstrack);
 910     nTrks = (Int_t)tracktree->GetEntries();
 911     nUsedTrks = 0;
 912     zvtxTrks = 0.;
 913     ezvtxTrks = 0.;
 914     sumWeights = 0.;
 915
 916     //cerr<<" ITS tracks: "<<nTrks<<endl;
 917     // loop on tracks
 918     for(i=0; i<nTrks; i++) {
 919       tracktree->GetEvent(i);
 920       // get pt and impact parameters
 921       pt = 1./TMath::Abs(itstrack->Get1Pt());
 922       d0rphi = TMath::Abs(itstrack->GetD());
 923       itstrack->PropagateToVertex();
 924       d0z = itstrack->GetZ();
 925       // check if the track is from (0,0) in (x,y)
 926       if(!VtxTrack(pt,d0rphi)) continue;
 927       nUsedTrks++;
 928       zvtxTrks   += d0z/d0zRes(pt)/d0zRes(pt);
 929       sumWeights += 1./d0zRes(pt)/d0zRes(pt);
 930
 931     } // loop on tracks
 932
 933     if(nUsedTrks>1) {
 934       // estimated position in z of vertex
 935       zvtxTrks /= sumWeights;
 936       // estimated error
 937       ezvtxTrks = TMath::Sqrt(1./sumWeights);
 938     } else {
 939       zvtxTrks = zvtxSPD[ev];
 940     }
 941
 942     zvtx[ev] = zvtxTrks;
 943
 944   } // loop on events
 945
 946
 947   // Write vertices to file
 948   WriteVtx("vtxTracks.root",zvtx,n);
 949   delete [] zvtx;
 950   delete [] zvtxSPD;
 951
 952   gBenchmark->Stop(name);
 953   gBenchmark->Show(name);
 954
 955   return 0;
 956 }
 957 //----------------------------------------------------------------------------
 958 Bool_t VtxTrack(Double_t pt,Double_t d0rphi) {
 959
 960   Double_t d0rphiRes = TMath::Sqrt(11.59*11.59+65.76*65.76/TMath::Power(pt,1.878))/10000.;
 961   Double_t nSigma = 3.;
 962   if(d0rphi<nSigma*d0rphiRes) { return kTRUE; } else { return kFALSE; }
 963 }
 964 //----------------------------------------------------------------------------
 965 Double_t d0zRes(Double_t pt) {
 966   Double_t res = TMath::Sqrt(34.05*34.05+170.1*170.1/TMath::Power(pt,1.226));
 967   return res/10000.;
 968 }
 969 //-----------------------------------------------------------------------------
 970 Int_t ITSFindTracksV2(const Char_t *galiceName,const Char_t *inName,
 971                       const Char_t *inName2,const Char_t *outName,
 972                       Int_t *skipEvt,TString vtxMode,Int_t n) {
 973
 974
 975   cerr<<"\n*******************************************************************\n";
 976
 977   const Char_t *name="ITSFindTracksV2";
 978   cerr<<'\n'<<name<<"...\n";
 979   gBenchmark->Start(name);
 980
 981   // Read vertices from file
 982   Double_t *zvtx = new Double_t[n];
 983   Char_t *vtxfile="vtxHeader.root";
 984
 985   cerr<<" Using vtxMode: ";
 986   switch (*(vtxMode.Data())) {
 987   case 'H':
 988     cerr<<" Header"<<endl;
 989     vtxfile = "vtxHeader.root";
 990     break;
 991   case 'F':
 992     cerr<<" fast pp"<<endl;
 993     vtxfile = "vtxFastpp.root";
 994     break;
 995   case 'P':
 996     cerr<<" SPD"<<endl;
 997     vtxfile = "vtxSPD.root";
 998     break;
 999   case 'T':
1000     cerr<<" Tracks"<<endl;
1001     vtxfile = "vtxTracks.root";
1002     break;
1003   }
1004
1005   ReadVtx(vtxfile,zvtx,n);
1006
1007
1008   TFile *outFile = TFile::Open(outName,"recreate");
1009   TFile *inFile  = TFile::Open(inName);
1010   TFile *inFile2 = TFile::Open(inName2);
1011
1012   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
1013   if(!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
1014
1015   Int_t flag1stPass,flag2ndPass;
1016
1017   // open logfile for done events
1018   FILE *logfile = fopen("itstracking.log","w");
1019
1020   AliITStrackerV2 tracker(geom);
1021   // loop on events
1022   for(Int_t ev=0; ev<n; ev++){
1023     // write to logfile of begun events
1024     fprintf(logfile,"%d\n",ev);
1025
1026     if(skipEvt[ev]) continue;
1027     cerr<<" --- Processing event "<<ev<<" ---"<<endl;
1028
1029     // pass event number to tracker
1030     tracker.SetEventNumber(ev);
1031
1032     // set position of primary vertex
1033     Double_t vtx[3];
1034     vtx[0]=0.; vtx[1]=0.; vtx[2]=zvtx[ev];
1035
1036     flag1stPass=1; // vtx constraint
1037     flag2ndPass=0; // no vtx constraint
1038
1039     // no vtx constraint if vertex not found
1040     if(vtx[2]<-999.) {
1041       flag1stPass=0;
1042       vtx[2]=0.;
1043     }
1044
1045     cerr<<"+++\n+++ Setting primary vertex z = "<<vtx[2]<<
1046       " cm for ITS tracking\n+++\n";
1047     tracker.SetVertex(vtx);
1048
1049     // setup vertex constraint in the two tracking passes
1050     Int_t flags[2];
1051     flags[0]=flag1stPass;
1052     tracker.SetupFirstPass(flags);
1053     flags[0]=flag2ndPass;
1054     tracker.SetupSecondPass(flags);
1055
1056     // find the tracks
1057     tracker.Clusters2Tracks(inFile,outFile);
1058
1059   } // loop on events
1060
1061   fprintf(logfile,"%d\n",n); //this means all evts are successfully completed
1062   fclose(logfile);
1063
1064   inFile->Close();
1065   inFile2->Close();
1066   outFile->Close();
1067
1068   delete [] zvtx;
1069
1070   gBenchmark->Stop(name);
1071   gBenchmark->Show(name);
1072
1073   return 0;
1074 }
1075 //-----------------------------------------------------------------------------
1076 Int_t ITSMakeRefFile(const Char_t *galice,const Char_t *inname,
1077                      const Char_t *outname,Int_t *skipEvt,Int_t n) {
1078
1079
1080   cerr<<"\n*******************************************************************\n";
1081
1082   Int_t rc=0;
1083   const Char_t *name="ITSMakeRefFile";
1084   cerr<<'\n'<<name<<"...\n";
1085   gBenchmark->Start(name);
1086
1087
1088   TFile *out = TFile::Open(outname,"recreate");
1089   TFile *trk = TFile::Open(inname);
1090   TFile *kin = TFile::Open(galice);
1091
1092
1093   // Get gAlice object from file
1094   if(!(gAlice=(AliRun*)kin->Get("gAlice"))) {
1095     cerr<<"gAlice has not been found on galice.root !\n";
1096     return 1;
1097   }
1098
1099   Int_t label;
1100   TParticle *Part;
1101   TParticle *Mum;
1102   RECTRACK rectrk;
1103
1104
1105   for(Int_t ev=0; ev<n; ev++){
1106     if(skipEvt[ev]) continue;
1107     cerr<<" --- Processing event "<<ev<<" ---"<<endl;
1108
1109     gAlice->GetEvent(ev);
1110
1111     trk->cd();
1112
1113     // Tree with ITS tracks
1114     char tname[100];
1115     sprintf(tname,"TreeT_ITS_%d",ev);
1116
1117     TTree *tracktree=(TTree*)trk->Get(tname);
1118     if(!tracktree) continue;
1119     AliITStrackV2 *itstrack=new AliITStrackV2;
1120     tracktree->SetBranchAddress("tracks",&itstrack);
1121     Int_t nentr=(Int_t)tracktree->GetEntries();
1122
1123     // Tree for true track parameters
1124     char ttname[100];
1125     sprintf(ttname,"Tree_Ref_%d",ev);
1126     TTree *reftree = new TTree(ttname,"Tree with true track params");
1127     reftree->Branch("rectracks",&rectrk,"lab/I:pdg:mumlab:mumpdg:Vx/F:Vy:Vz:Px:Py:Pz");
1128
1129     for(Int_t i=0; i<nentr; i++) {
1130       tracktree->GetEvent(i);
1131       label = TMath::Abs(itstrack->GetLabel());
1132
1133       Part = (TParticle*)gAlice->Particle(label);
1134       rectrk.lab=label;
1135       rectrk.pdg=Part->GetPdgCode();
1136       rectrk.mumlab = Part->GetFirstMother();
1137       if(Part->GetFirstMother()>=0) {
1138         Mum = (TParticle*)gAlice->Particle(Part->GetFirstMother());
1139         rectrk.mumpdg=Mum->GetPdgCode();
1140       } else {
1141         rectrk.mumpdg=-1;
1142       }
1143       rectrk.Vx=Part->Vx();
1144       rectrk.Vy=Part->Vy();
1145       rectrk.Vz=Part->Vz();
1146       rectrk.Px=Part->Px();
1147       rectrk.Py=Part->Py();
1148       rectrk.Pz=Part->Pz();
1149
1150       reftree->Fill();
1151     } // loop on tracks
1152
1153     out->cd();
1154     reftree->Write();
1155
1156     delete itstrack;
1157     delete reftree;
1158   } // loop on events
1159
1160   trk->Close();
1161   kin->Close();
1162   out->Close();
1163
1164   gBenchmark->Stop(name);
1165   gBenchmark->Show(name);
1166
1167
1168   return rc;
1169 }
1170 //-----------------------------------------------------------------------------
1171 void WriteVtx(const Char_t *name,Double_t *zvtx,Int_t n) {
1172
1173   // Set Random Number seed
1174   TDatime dt;
1175   UInt_t curtime=dt.Get();
1176   UInt_t procid=gSystem->GetPid();
1177   UInt_t seed=curtime-procid;
1178   gRandom->SetSeed(seed);
1179
1180   Double_t xvtx,yvtx;
1181   Double_t sigmaVtxTransverse = 15.e-4;
1182   TVector3 *ioVtx = new TVector3;
1183   TTree *vtxtree = new TTree("TreeVtx","Tree with positions of primary vertex");
1184   vtxtree->Branch("vertices","TVector3",&ioVtx);
1185
1186   for(Int_t ev=0; ev<n; ev++) {
1187     // x and y coordinates of the vertex are smeared according to a gaussian
1188     xvtx = gRandom->Gaus(0.,sigmaVtxTransverse);
1189     yvtx = gRandom->Gaus(0.,sigmaVtxTransverse);
1190
1191     ioVtx->SetX(xvtx);
1192     ioVtx->SetY(yvtx);
1193     ioVtx->SetZ(zvtx[ev]);
1194
1195     //cerr<<"W ev "<<ev<<"  ("<<ioVtx->X()<<","<<ioVtx->Y()<<","<<ioVtx->Z()<<")"<<endl;
1196
1197
1198     vtxtree->Fill();
1199   }
1200
1201   // write the tree to a file
1202   TFile *f = new TFile(name,"recreate");
1203   vtxtree->Write();
1204   f->Close();
1205
1206   return;
1207 }
1208 //-----------------------------------------------------------------------------
1209 void ReadVtx(const Char_t *name,Double_t *zvtx,Int_t n) {
1210
1211   TFile *f = TFile::Open(name);
1212
1213   TVector3 *ioVtx = 0;
1214   TTree *vtxtree = (TTree*)f->Get("TreeVtx");
1215   vtxtree->SetBranchAddress("vertices",&ioVtx);
1216   Int_t entries = (Int_t)vtxtree->GetEntries();
1217
1218   for(Int_t ev=0; ev<n; ev++) {
1219     vtxtree->GetEvent(ev);
1220
1221     zvtx[ev] = ioVtx->Z();
1222
1223     //cerr<<"R ev "<<ev<<"  ("<<ioVtx->X()<<","<<ioVtx->Y()<<","<<zvtx[ev]<<")"<<endl;
1224
1225   }
1226
1227   f->Close();
1228
1229   return;
1230 }
1231 //-----------------------------------------------------------------------------
1232
1233
1234
1235
1236
1237
1238
1239
1240