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