ITS/AliITSandTPCrecoV2.C

   1 ////////////////////////////////////////////////////////////////////////
   2 //
   3 // AliITSandTPCrecoV2.C
   4 //
   5 //
   6 ////////////////////////////////////////////////////////////////////////
   7
   8
   9 /****************************************************************************
  10  * This macro is supposed to do reconstruction in the barrel ALICE trackers *
  11  * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
  12  * It does the following steps                                              *
  13  *                   1) TPC cluster finding                                 *
  14  *                   2) TPC track finding                                   *
  15  *                   3) ITS cluster finding - fast recpoints -              *
  16  *                   4) ITS tracking                                        *
  17  * TPC part taken from AliTPCTracking.C - Jiri Chudoba                      *
  18  * The TPC part will be removed                                             *
  19  * M.Masera 30/09/2002 09:40                                                *
  20  * Last revision: 4/3/2003                                                  *
  21  *                                                                          *
  22  * Use:                                                                     *
  23  * The main function is AliITSandTPCrecoV2 All the input arguments have a   *
  24  * default value.                                                           *
  25  * nEvents: number of events to be processed.                               *
  26  *          If nEvents<0, all the available events are used starting from   *
  27  *          the first one                                                   *
  28  * firstEvent: is the first event to be analyzed.                           *
  29  *             It is set to 0 if nEvents<0 (default)                        *
  30  * fileNameHits: input root file with Kine and hits                         *
  31  * fileNameDigits: input root file with ITS and TPC digits and with ITS     *
  32  *                 recpoints. Only fast points are used at the moment.      *
  33  *                 An option to use slow points should be added             *
  34  * fileNameClusters: output root file with TPC clusters                     *
  35  * fileNameTracks: output root file with TPC tracks                         *
  36  * fileNameITSClusters: output root file with ITS clusters                  *
  37  * fileNameITSTracks: output root file with ITS tracks                      *
  38  *                                                                          *
  39  *   AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0,               *
  40  *                   TString fileNameHits="galice.root",                    *
  41  *                   TString fileNameDigits="galiceSDR.root",               *
  42  *                   TString fileNameClusters="tpc.clusters.root",          *
  43  *                   TString fileNameTracks="tpc.tracks.root",              *
  44  *                   TString fileNameITSClusters="its.clusters.root",       *
  45  *                   TString fileNameITSTracks="its.tracks.root");          *
  46  *                                                                          *
  47  *  Global variables:
  48  *
  49  * Int_t gDEBUG = 0;   // NO verbose printouts if ==0                       *
  50  * TArrayD *gzvtx = 0; // Z coordinate of primary vertices                  *
  51  * Int_t gCurrEve = 0; // current event number                              *
  52  * Bool_t gPPMode = kTRUE;  // it must be set to kTRUE for pp interactions  *
  53  * Bool_t gFastReco = kTRUE;  // it must set to kTRUE for fast reconstr.    *
  54  * const Double_t kgVertexNotFound = -100;  // the Z coord. of the          *
  55  *                                 primary vertex is set to this value      *
  56  *                                 when the vertexing fails                 *
  57  * Bool_t gUseNewClustersV2 = kFALSE;  // if kTRUE the AliITSclustereV2 class
  58  *                                 is used to produce V2 clusters.          *
  59  *                                 ITS recpoints are used otherwise         *
  60  ****************************************************************************/
  61
  62 #if !defined(__CINT__) || defined(__MAKECINT__)
  63 #include <TH1.h>
  64 #include <TH2.h>
  65 #include <TProfile.h>
  66 #include <TBranch.h>
  67 #include <TParticle.h>
  68 #include <TRandom.h>
  69 #include "alles.h"
  70 #include "AliMagF.h"
  71 #include "AliMagFCM.h"
  72 #include "AliTPCtracker.h"
  73 #include "AliTPCclusterer.h"
  74 #include "Riostream.h"
  75 #include "TString.h"
  76 #include "TArrayD.h"
  77 #include "TLine.h"
  78 #include "TTree.h"
  79 #include "TNtuple.h"
  80 #include "TParticle.h"
  81 #include "TText.h"
  82 #include "TSystem.h"
  83 #include "TVector.h"
  84 #include "TObjArray.h"
  85 #include "TStyle.h"
  86 #include "TCanvas.h"
  87 #include "AliITS.h"
  88 #include "AliITSclustererV2.h"
  89 #include "AliITSgeom.h"
  90 #include "AliITSRecPoint.h"
  91 #include "AliITSclusterV2.h"
  92 #include "AliITStrackerV2.h"
  93 #include "AliITSVertexerPPZ.h"
  94 #include "AliITSVertexerIons.h"
  95 #include "AliRun.h"
  96 #include "AliHeader.h"
  97 #include "AliGenEventHeader.h"
  98 #include "TArrayF.h"
  99 #include "AliGenPythiaEventHeader.h"
 100
 101 #endif
 102
 103
 104 Int_t gDEBUG = 0;
 105 TArrayD *gzvtx = 0;
 106 Int_t gCurrEve = 0;
 107 Bool_t gPPMode = kTRUE;
 108 Bool_t gFastReco = kTRUE;
 109 Bool_t gUseNewClustersV2 = kFALSE;
 110 const Double_t kgVertexNotFound = -100.;
 111
 112
 113 Int_t TPCFindClusters(Int_t nEvents=1, Int_t firstEvent=0,
 114                       TString fileNameDigits="galiceSDR.root",
 115                       TString fileNameClusters="tpc.clusters.root");
 116 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
 117                       TFile* fileDigits, TFile* fileClusters,
 118                       AliTPCParam* paramTPC=0);
 119 Int_t TPCFindTracks(Int_t nEvents=1, Int_t firstEvent=0,
 120                     TString fileNameClusters="tpc.clusters.root",
 121                     TString fileNameTracks="tpc.tracks.root");
 122 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
 123                     TFile* fileClusters, TFile* fileTracks,
 124                     AliTPCParam* paramTPC=0);
 125
 126 Int_t AliTPCTracking(Int_t nEvents=1, Int_t firstEvent=0,
 127                      TString fileNameHits="galice.root",
 128                      TString fileNameDigits="galiceSDR.root",
 129                      TString fileNameClusters="tpc.clusters.root",
 130                      TString fileNameTracks="tpc.tracks.root");
 131 Int_t AliITSandTPCrecoV2(Int_t nEvents=-1, Int_t firstEvent=0,
 132                         TString fileNameHits="galice.root",
 133                         TString fileNameDigits="galiceSDR.root",
 134                         TString fileNameClusters="tpc.clusters.root",
 135                         TString fileNameTracks="tpc.tracks.root",
 136                         TString fileNameITSClusters="its.clusters.root",
 137                         TString fileNameITSTracks="its.tracks.root");
 138 Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first);
 139 Int_t ITSFindTracks(TString fileNameTracks, TString fileNameITSClusters, TString fileNameITSTracks, Int_t n, Int_t first, Int_t vc, Int_t vc2);
 140 void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints);
 141 Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile *in,TFile *out);
 142 ////////////////////////////////////////////////////////////////////////
 143 Int_t AliITSandTPCrecoV2(Int_t nEvents, Int_t firstEvent,
 144                         TString fileNameHits,
 145                         TString fileNameDigits,
 146                         TString fileNameClusters,
 147                         TString fileNameTracks,
 148                         TString fileNameITSClusters,
 149                         TString fileNameITSTracks) {
 150   const Char_t *name="AliITSandTPCrecoV2";
 151
 152   gBenchmark->Start(name);
 153
 154   // Set the conversion constant for the Kalman filter
 155   // and set the gAlice pointer
 156   Char_t *finame = (Char_t *)fileNameHits.Data();
 157   AliTracker::SetFieldFactor(finame,kFALSE);
 158
 159   if(nEvents < 0 ) {
 160     nEvents = (Int_t)gAlice->TreeE()->GetEntries();
 161     firstEvent = 0;
 162     cout << "Processing events from " << firstEvent << " up to " << nEvents-1 <<endl;
 163   }
 164
 165   if(fileNameDigits!=fileNameHits)gAlice->SetTreeDFileName(fileNameDigits.Data());
 166   gzvtx = new TArrayD(nEvents);
 167
 168   // measure Z vertex coordinate for seeding
 169
 170   FindZV(nEvents,firstEvent,fileNameHits,fileNameDigits);
 171
 172   // TPC tracking - abort macro execution if tracking fails
 173   if(AliTPCTracking(nEvents,firstEvent,fileNameHits,fileNameDigits,
 174                     fileNameClusters,fileNameTracks)) {
 175     cerr << "TPC tracking failed \n";
 176     return 1;
 177   }
 178
 179   // ITS clustering
 180   if(ITSFindClusters(fileNameHits,fileNameDigits,fileNameITSClusters,nEvents,firstEvent)) {
 181     cerr << "ITS clustering failed \n";
 182     return 2;
 183   }
 184
 185   Int_t firstpass = 1;  // 1=vert. constraint; 0=no vert. constr.; -1=skip pass
 186   Int_t secondpass = 0; // as above
 187   if(ITSFindTracks(fileNameTracks,fileNameITSClusters,fileNameITSTracks,nEvents,firstEvent,firstpass,secondpass)) {
 188     cerr << "ITS tracking failed \n";
 189     return 3;
 190   }
 191
 192
 193   gBenchmark->Stop(name);
 194   gBenchmark->Show(name);
 195   return 0;
 196 }
 197 ////////////////////////////////////////////////////////////////////////
 198 Int_t AliTPCTracking( Int_t nEvents, Int_t firstEvent,
 199                       TString fileNameHits,
 200                       TString fileNameDigits,
 201                       TString fileNameClusters,
 202                       TString fileNameTracks) {
 203
 204
 205   // ********** Find TPC clusters *********** //
 206   if (TPCFindClusters(nEvents,firstEvent,fileNameDigits,fileNameClusters)) {
 207     cerr<<"Failed to get TPC clusters: !\n";
 208     return 1;
 209   }
 210
 211   // ********** Find TPC tracks *********** //
 212   //   if (TPCFindTracks(TPCclsName,TPCtrkName,n)) {
 213   if (TPCFindTracks(nEvents,firstEvent,fileNameClusters,fileNameTracks)) {
 214     cerr<<"Failed to get TPC tracks !\n";
 215     return 2;
 216   }
 217
 218   return 0;
 219 }
 220
 221 ////////////////////////////////////////////////////////////////////////
 222 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
 223                       TString fileNameDigits, TString fileNameClusters) {
 224
 225   Int_t rc;
 226   const Char_t *name="TPCFindClusters";
 227   if (gDEBUG>1) cout<<name<<" starts...\n";
 228   if (gDEBUG>1) gBenchmark->Start(name);
 229   TFile *fileClusters = TFile::Open(fileNameClusters,"recreate");
 230   TFile *fileDigits = TFile::Open(fileNameDigits);
 231   if (!fileDigits->IsOpen()) {
 232     cerr<<"Cannnot open "<<fileNameDigits<<" !\n";
 233     return 1;
 234   }
 235   if (!fileClusters->IsOpen()) {
 236     cerr<<"Cannnot open "<<fileNameClusters<<" !\n";
 237     return 1;
 238   }
 239
 240   rc = TPCFindClusters(nEvents,firstEvent,fileDigits,fileClusters);
 241
 242   fileDigits->Close();
 243   fileClusters->Close();
 244   delete fileDigits;
 245   delete fileClusters;
 246   if (gDEBUG>1) gBenchmark->Stop(name);
 247   if (gDEBUG>1) gBenchmark->Show(name);
 248
 249   return rc;
 250 }
 251 ////////////////////////////////////////////////////////////////////////
 252 Int_t TPCFindClusters(Int_t nEvents, Int_t firstEvent,
 253                       TFile* fileDigits, TFile* fileClusters,
 254                       AliTPCParam* paramTPC) {
 255
 256   fileDigits->cd();
 257   if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileDigits);
 258   if (!paramTPC) return 1;
 259
 260   for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
 261     if (gDEBUG > 2) cout<<"TPCFindClusters: event "<<iEvent<<endl;
 262     AliTPCclusterer::Digits2Clusters(paramTPC, fileClusters, iEvent);
 263   }
 264   return 0;
 265 }
 266 ////////////////////////////////////////////////////////////////////////
 267 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
 268                     TString fileNameClusters, TString fileNameTracks) {
 269
 270   Int_t rc = 0;
 271   const Char_t *name="TPCFindTracks";
 272   if (gDEBUG>1) cout<<name<<" starts"<<endl;
 273   if (gDEBUG>1) gBenchmark->Start(name);
 274   TFile *fileTracks = TFile::Open(fileNameTracks,"recreate");
 275   TFile *fileClusters =TFile::Open(fileNameClusters);
 276
 277   rc = TPCFindTracks(nEvents, firstEvent, fileClusters, fileTracks);
 278
 279   fileClusters->Close();
 280   fileTracks->Close();
 281   delete fileClusters;
 282   delete fileTracks;
 283   if (gDEBUG>1) gBenchmark->Stop(name);
 284   if (gDEBUG>1) gBenchmark->Show(name);
 285   return rc;
 286
 287 }
 288 ////////////////////////////////////////////////////////////////////////
 289 Int_t TPCFindTracks(Int_t nEvents, Int_t firstEvent,
 290                     TFile *fileClusters, TFile * fileTracks,
 291                     AliTPCParam* paramTPC) {
 292
 293   Int_t rc = 0;
 294   if (!paramTPC) paramTPC = AliTPC::LoadTPCParam(fileClusters);
 295   if (!paramTPC) return 1;
 296   Double_t vertex[3];
 297   for(Int_t j=0; j<3; j++){vertex[j]=0.;}
 298   AliTPCtracker *tracker = new AliTPCtracker(paramTPC);
 299   for (Int_t iEvent = firstEvent; iEvent < firstEvent+nEvents; iEvent++){
 300     if (gDEBUG > 2) cout<<"TPCFindTracks: event "<<iEvent<<endl;
 301     tracker->SetEventNumber(iEvent);
 302     if((*gzvtx)[iEvent] > -100){
 303       vertex[2] = (*gzvtx)[iEvent];
 304     }
 305     else {
 306       vertex[2]=0;
 307     }
 308     tracker->SetVertex(vertex);
 309     fileClusters->cd();
 310     rc = tracker->Clusters2Tracks(0,fileTracks);
 311   }
 312   delete tracker;
 313   return rc;
 314 }
 315
 316
 317 ////////////////////////////////////////////////////////////////////////
 318 Int_t ITSFindClusters(TString inname, TString recpointFileName, TString fileNameITSClusters, Int_t n, Int_t first) {
 319   Int_t rc=0;
 320   const Char_t *name="ITSFindClusters";
 321   cerr<<'\n'<<name<<"...\n";
 322   gBenchmark->Start(name);
 323   TFile *out=TFile::Open(fileNameITSClusters,"recreate");
 324   TFile *in = (TFile*)gROOT->GetListOfFiles()->FindObject(inname.Data());
 325   if(gUseNewClustersV2){
 326     cout<<"Direct method\n";
 327     return DirectClusterFinder(n,first,in,out);
 328   }
 329   else {
 330     cout<<"Clusters through rec points \n";
 331   }
 332
 333   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 334   if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
 335   AliITSgeom *geom=ITS->GetITSgeom();
 336   AliITSclustererV2 *clusterer = new AliITSclustererV2(geom);
 337   out->cd();
 338   geom->Write();
 339
 340   Int_t ev;
 341   if(gFastReco && gDEBUG>0)cout <<"Using fast points\n";
 342   else cout<<"Using slow points\n";
 343   for (ev = first; ev<n; ev++){
 344     cout<<"ITSFindClusters: processing event "<<ev<<endl;
 345     in->cd();
 346     gAlice->GetEvent(ev);
 347
 348     TTree *pTree=gAlice->TreeR();
 349     if (!pTree){
 350       cerr << "ITSFindClusters: can't get TreeR\n";
 351       return 1;
 352     }
 353     TBranch *branch = 0;
 354     if (gFastReco) {
 355       branch = pTree->GetBranch("ITSRecPointsF");
 356     }
 357     else {
 358       branch = pTree->GetBranch("ITSRecPoints");
 359     }
 360     if (!branch) {
 361       cerr << "ITSFindClusters: can't get ITSRecPoints branch \n";
 362       return 2;
 363     }
 364
 365     out->cd();
 366     TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
 367     char   cname[100];
 368     sprintf(cname,"TreeC_ITS_%d",ev);
 369
 370     TTree *cTree=new TTree(cname,"ITS clusters");
 371     cTree->Branch("Clusters",&clusters);
 372
 373     TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
 374     branch->SetAddress(&points);
 375
 376     TClonesArray &cl=*clusters;
 377     Int_t nclusters=0;
 378     Int_t nentr=(Int_t)pTree->GetEntries();
 379     AliITSgeom *geom=ITS->GetITSgeom();
 380     Float_t lp[5];
 381     Double_t rot[9];
 382     Int_t lab[6];
 383     for (Int_t i=0; i<nentr; i++) {
 384       branch->GetEvent(i);
 385       clusterer->RecPoints2Clusters(points,i,clusters);
 386       cTree->Fill(); clusters->Delete();
 387       points->Clear();
 388     }
 389     cTree->Write();
 390     delete cTree; delete clusters; points->Delete(); delete points;
 391   }
 392   delete clusterer;
 393
 394   out->Close();
 395
 396   return rc;
 397 }
 398 ////////////////////////////////////////////////////////////////////////
 399 Int_t ITSFindTracks(TString inname, TString inname2, TString outname, Int_t n, Int_t first, Int_t vc, Int_t vc2) {
 400   Int_t rc=0;
 401   if(gPPMode){
 402     AliITStrackV2::SetSigmaYV(0.015);
 403     AliITStrackV2::SetSigmaZV(0.017);
 404   }
 405   TFile *out=TFile::Open(outname.Data(),"recreate");
 406   if (!out->IsOpen()) {
 407     cerr<<"Can't open file "<<outname.Data()<<endl;
 408     return 1;
 409   }
 410   TFile *in =TFile::Open(inname.Data());
 411   if (!in->IsOpen()) {
 412     cerr<<"Can't open file "<<inname.Data()<<endl;
 413     return 1;
 414   }
 415   TFile *in2 =TFile::Open(inname2.Data());
 416   if (!in2->IsOpen()) {
 417     cerr<<"Can't open file "<<inname2.Data()<<endl;
 418     return 1;
 419   }
 420
 421   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
 422   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
 423   Double_t vrtc[3];
 424   AliITStrackerV2 tracker(geom);
 425   for (Int_t i=first;i<n;i++){
 426     tracker.SetEventNumber(i);
 427     Double_t vtx=(*gzvtx)[i];
 428     cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++\n";
 429     cout<<"ITSFindTracks -- event "<<i<<endl;
 430     cout<<"Computed vertex position: "<<vtx<<endl;
 431     in2->cd();
 432     vrtc[0]=0.;
 433     vrtc[1]=0.;
 434     vrtc[2]=vtx;
 435     if(vtx != kgVertexNotFound){
 436       tracker.SetVertex(vrtc);   // set vertex only if it was computed
 437     }
 438     else {
 439       vc = 0;  // use vertex contraint only if vertex info is available
 440       vc2 = -1;
 441     }
 442
 443     tracker.SetupFirstPass(&vc);
 444
 445     tracker.SetupSecondPass(&vc2);
 446
 447     rc=tracker.Clusters2Tracks(in,out);
 448   }
 449
 450   in->Close();
 451   in2->Close();
 452   out->Close();
 453
 454   return rc;
 455 }
 456
 457 //////////////////////////////////////////////////////////////////
 458 void FindZV(Int_t nEvents, Int_t first, TString FileNameHits, TString FileWithRecPoints){
 459   TFile *in = (TFile*)gROOT->GetListOfFiles()->FindObject(FileNameHits.Data());
 460   TFile *fo = new TFile("AliITSVertices.root","recreate");
 461   if(FileNameHits!=FileWithRecPoints)gAlice->SetTreeRFileName(FileWithRecPoints);
 462   AliITSVertexer *findVert;
 463   if(gPPMode){
 464     findVert = new AliITSVertexerPPZ(in,fo);
 465   }
 466   else {
 467     findVert = new AliITSVertexerIons(in,fo);
 468   }
 469   Int_t last = first + nEvents - 1;
 470   AliITSVertex *vert = 0;
 471   for(Int_t i=first; i<=last; i++){
 472     gAlice->GetEvent(i);
 473     // The true Z coord. is fetched for comparison
 474     AliHeader *header = gAlice->GetHeader();
 475     AliGenEventHeader* genEventHeader = header->GenEventHeader();
 476     TArrayF primaryVertex(3);
 477     genEventHeader->PrimaryVertex(primaryVertex);
 478     vert = findVert->FindVertexForCurrentEvent(i);
 479     if(vert){
 480       Double_t pos[3];
 481       for(Int_t kk=0;kk<3;kk++)pos[kk]=(Double_t)primaryVertex[kk];
 482       vert->SetTruePos(pos);
 483       Double_t found = vert->GetZv();
 484       gzvtx->AddAt(found,i);
 485       findVert->WriteCurrentVertex();
 486     }
 487     else {
 488       gzvtx->AddAt(kgVertexNotFound,i);
 489     }
 490   }
 491   delete findVert;
 492   fo->Close();
 493   delete fo;
 494 }
 495
 496
 497
 498 Int_t DirectClusterFinder(Int_t eventn,Int_t first,TFile *in,TFile* out){
 499
 500
 501    AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
 502    if (!ITS) { cerr<<"Can't find the ITS !\n"; return 3; }
 503
 504    AliITSgeom *geom=ITS->GetITSgeom();
 505
 506    geom->Write();
 507
 508    TStopwatch timer;
 509    AliITSclustererV2 clusterer(geom);
 510    for (Int_t i=first; i<eventn; i++) {
 511        cerr<<"Processing event number: "<<i<<endl;
 512        gAlice->GetEvent(i);
 513        //ITS->MakeTreeC(); //To make the V1 cluster finders happy
 514        clusterer.SetEvent(i);
 515        if (gFastReco) clusterer.Digits2Clusters(in,out);
 516        else                 clusterer.Hits2Clusters(in,out);
 517    }
 518    timer.Stop(); timer.Print();
 519
 520    out->Close();
 521
 522    return 0;
 523
 524 }
 525