X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliBarrelRec_TPCparam.C;h=5bf2da2118594b1a012cc9856c7b2311f6aabddd;hb=3f3beb83dfda9261f00e33673d047d5acaede4db;hp=b04dd08c3dc3e64f98d8911b3dcc53484003154a;hpb=a654e2576b86596dc5a0d70b90fdc22849fd3ac3;p=u%2Fmrichter%2FAliRoot.git

diff --git a/TPC/AliBarrelRec_TPCparam.C b/TPC/AliBarrelRec_TPCparam.C
index b04dd08c3dc..5bf2da21185 100644
--- a/TPC/AliBarrelRec_TPCparam.C
+++ b/TPC/AliBarrelRec_TPCparam.C
@@ -1,7 +1,7 @@
 void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
   //
   // Macro to create AliESDs.root using parametrized TPC tracking
-  // and AliITStrackerV2
+  // and AliITStrackerSA (MI + SA)
   //
   // Input files:
   // - galice.root
@@ -10,155 +10,244 @@ void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
   // - ITS.RecPoints.root (use AliRecontruction class)
   // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
   //
-  // A. Dainese - LNL
+  // A. Dainese - INFN Legnaro
   //
 
 
-   /**** Initialization of the NewIO *******/
-
-   if (gAlice) {
-      delete gAlice->GetRunLoader();
-      delete gAlice; 
-      gAlice=0;
-   }
-
-   AliRunLoader *rl = AliRunLoader::Open("galice.root");
-   if (rl == 0x0) {
-      cerr<<"Can not open session"<<endl;
-      return;
-   }
-   Int_t retval = rl->LoadgAlice();
-   if (retval) {
-      cerr<<"LoadgAlice returned error"<<endl;
-      delete rl;
-      return;
-   }
-   retval = rl->LoadHeader();
-   if (retval) {
-      cerr<<"LoadHeader returned error"<<endl;
-      delete rl;
+  Int_t  collcode = 1; // pp collisions
+  Bool_t useMeanVtx = kFALSE;
+  
+  AliGeomManager::LoadGeometry("geometry.root");
+  
+  if (gAlice) {
+    delete gAlice->GetRunLoader();
+    delete gAlice; 
+    gAlice=0;
+  }
+  AliRunLoader *rl = AliRunLoader::Open("galice.root");
+  if (rl == 0x0) {
+    cerr<<"Can not open session"<<endl;
+    return;
+  }
+  Int_t retval = rl->LoadgAlice();
+  if (retval) {
+    cerr<<"LoadgAlice returned error"<<endl;
+    delete rl;
+    return;
+  }
+  retval = rl->LoadHeader();
+  if (retval) {
+    cerr<<"LoadHeader returned error"<<endl;
+    delete rl;
+    return;
+  }
+  gAlice=rl->GetAliRun();
+  
+  
+  
+  // Get field from galice.root
+  AliMagF *fiel = (AliMagF*)gAlice->Field();
+  Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
+  // Set the conversion constant between curvature and Pt
+  AliTracker::SetFieldMap(fiel,kTRUE);
+  
+  /**** The TPC corner ********************/
+  
+  AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
+  tpcTrackerPar.Init();
+  
+  //**** Switch on the PID class (mandatory!)
+  AliPID pid;
+  
+  /**** The ITS corner ********************/
+  
+  AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
+  if (itsl == 0x0) {
+    cerr<<"Cannot get the ITS loader"<<endl;
+    return;
+  }
+  itsl->LoadRecPoints("read");
+  
+  AliITSRecoParam * itsRecoParam = AliITSRecoParam::GetLowFluxParam();
+  AliITSReconstructor::SetRecoParam(itsRecoParam);
+  
+  // Instance of the ITS tracker
+  AliITStrackerSA itsTracker(0);
+  Int_t ITSclusters[6] = {1,1,1,1,1,1};
+  itsTracker.SetLayersNotToSkip(ITSclusters);  
+  
+  // Primary vertex reconstruction in pp
+  AliESDVertex *initVertex = 0;
+  TFile *invtx = new TFile("AliESDVertexMean.root");
+  if(collcode==1 && useMeanVtx) {
+    // open the mean vertex
+    initVertex = (AliESDVertex*)invtx->Get("vtxmean");
+  } else {
+    Double_t pos[3]={0.5,0.5,0.}; 
+    Double_t err[3]={3.,3.,5.};
+    initVertex = new AliESDVertex(pos,err);
+  }
+  invtx->Close();
+  delete invtx;
+  AliVertexerTracks *vertexer = new AliVertexerTracks(AliTracker::GetBz());
+  vertexer->SetVtxStart(initVertex);
+  vertexer->SetDebug(0);
+  delete initVertex;
+  initVertex=0;
+  
+  /***** The TREE for ESD is born *****/
+  TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
+  AliESDEvent *event=0; AliESDEvent *eventTPCin=0;
+  event = new AliESDEvent();
+  event->CreateStdContent();
+  event->WriteToTree(esdTree);
+  
+  if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
+  if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
+  cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
+  
+  TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
+  AliESDVertex *vertexSPD = new AliESDVertex();
+  Char_t zver[100];
+  Double_t vtx[3]={0,0,0};
+  Double_t sigmavtx[3]={0.07,0.07,0.1};
+
+
+  //<---------------------------------- The Loop over events begins
+  TStopwatch timer;
+  Int_t trc;
+  for(Int_t i=firstEvent; i<=lastEvent; i++) { 
+    
+    cout<<" Processing event number : "<<i<<endl;
+    //AliESDEvent *event = new AliESDEvent(); 
+    event->SetRunNumber(gAlice->GetRunNumber());
+    event->SetEventNumberInFile(i);
+    event->SetMagneticField(gAlice->Field()->SolenoidField());
+    rl->GetEvent(i);
+
+    //***** Primary vertex from SPD from file 
+    sprintf(zver,"Event%d/Vertex",i);
+    vertexSPD = (AliESDVertex*)ppZ->Get(zver);
+    if(!vertexSPD) {
+      esdTree->Fill(); event->Reset();
+      continue;
+    }      
+    event->SetVertex(vertexSPD);
+    vertexSPD->GetXYZ(vtx);
+    vertexSPD->GetSigmaXYZ(sigmavtx);
+
+    //***** TPC tracking
+    if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
+      printf("exiting TPC tracker with code %d in event %d\n",trc,i);
+      esdTree->Fill(); event->Reset();
+      continue;
+    }
+
+    // make a copy of the ESD at this stage
+    eventTPCin = event;
+
+    //***** ITS tracking
+    itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
+    //    itsl->LoadRecPoints("read");
+    TTree *itsTree=itsl->TreeR();
+    if (!itsTree) {
+      cerr<<"Can't get the ITS cluster tree !\n";
+      esdTree->Fill(); event->Reset();
       return;
-   }
-   gAlice=rl->GetAliRun();
-       
-
-   TDatabasePDG *DataBase = TDatabasePDG::Instance();
-
-   // Get field from galice.root
-   AliMagF *fiel = (AliMagF*)gAlice->Field();
-   Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
-   // Set the conversion constant between curvature and Pt
-   AliTracker::SetFieldMap(fiel,kTRUE);
-
-   /**** The TPC corner ********************/
-
-   Int_t collcode = 1; // pp collisions
-   AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
-   tpcTrackerPar.Init();
-
-   //**** Switch on the PID class (mandatory!)
-   AliPID pid;
-
-   /**** The ITS corner ********************/
-
-   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
-   if (itsl == 0x0) {
-      cerr<<"Cannot get the ITS loader"<<endl;
-      return;
-   }
-   itsl->LoadRecPoints("read");
-
-   AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
-   if (!dITS) {
-      cerr<<"Cannot find the ITS detector !"<<endl;
-      return;
-   }
-   AliITSgeom *geom = dITS->GetITSgeom();
-
-   //An instance of the ITS tracker
-   AliITStrackerV2 itsTracker(geom);
-   Int_t ITSclusters[6] = {1,1,1,1,1,1};
-   itsTracker.SetLayersNotToSkip(ITSclusters);
-
-   /***** The TREE is born *****/
-   
-   TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
-   AliESD *event=0;
-   esdTree->Branch("ESD","AliESD",&event);
-   
-   if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
-   if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
-   cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
-   
-   TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
-   AliESDVertex *myvertex = new AliESDVertex();
-   Char_t zver[100];
-   Double_t vtx[3];
-   Double_t cvtx[6];
-
-
-   //<----------------------------------The Loop over events begins
-   TStopwatch timer;
-   Int_t trc;
-   for(Int_t i=firstEvent; i<=lastEvent; i++) { 
-     
-     cerr<<" Processing event number : "<<i<<endl;
-     AliESD *event = new AliESD(); 
-     event->SetRunNumber(gAlice->GetRunNumber());
-     event->SetEventNumber(i);
-     event->SetMagneticField(gAlice->Field()->SolenoidField());
-     rl->GetEvent(i);
-
-    //***** Primary vertex 
-     sprintf(zver,"Event%d/Vertex",i);
-     myvertex = (AliESDVertex*)ppZ->Get(zver);
-     if(!myvertex) {
-       esdTree->Fill(); delete event;
-       continue;
-     }      
-     event->SetVertex(myvertex);
-     myvertex->GetXYZ(vtx);
-     myvertex->GetCovMatrix(cvtx);
-
-     //***** TPC tracking
-     if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
-       printf("exiting TPC tracker with code %d in event %d\n",trc,i);
-       esdTree->Fill(); delete event;
-       continue;
-     }
-
-     //***** ITS tracking
-     itsTracker.SetVertex(vtx,cvtx);
-     TTree *itsTree=itsl->TreeR();
-     if (!itsTree) {
-        cerr<<"Can't get the ITS cluster tree !\n";
-	esdTree->Fill(); delete event;
-        return;
-     }     
-     itsTracker.UnloadClusters();
-     itsTracker.LoadClusters(itsTree);
-     if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
-       printf("exiting ITS tracker with code %d in event %d\n",trc,i);
-       esdTree->Fill(); delete event;
-       continue;
-     }
-
-
-     esdTree->Fill();
-     delete event;
-
-   }//<-----------------------------------The Loop over events ends here
-   timer.Stop(); timer.Print();
-
-   //        The AliESDs.root is born
-   TFile *ef = TFile::Open("AliESDs.root","RECREATE"); 
-   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
-
-   esdTree->Write(); //Write the TREE and close everything
-   delete esdTree;
-   ef->Close();
-
-   delete rl;
-
-   return;
+    }     
+    itsTracker.UnloadClusters();
+    itsTracker.LoadClusters(itsTree);
+    if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
+      printf("exiting ITS tracker with code %d in event %d\n",trc,i);
+      esdTree->Fill(); event->Reset();
+      continue;
+    }
+
+    // Bring kTPCin-tracks back to the TPC inner wall
+    BackToTPCInnerWall(event,eventTPCin);
+
+    // refit inward in ITS:
+    // - refit without vertex constraint
+    // - propagate through beam pipe to local x = 0
+    itsTracker.RefitInward(event);
+    
+
+    //***** Vertex from ESD tracks
+    if(collcode==1) { // pp
+      AliESDVertex *vertexTrks = 
+	(AliESDVertex*)vertexer->FindPrimaryVertex(event);
+      event->SetPrimaryVertex(vertexTrks); 
+    }
+    
+    esdTree->Fill();
+    event->Reset();
+
+  }//<-----------------------------------The Loop over events ends here
+  timer.Stop(); timer.Print();
+  
+  // The AliESDs.root is born
+  TFile *ef = TFile::Open("AliESDs.root","RECREATE"); 
+  if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
+
+  //Write the tree and close everything
+  esdTree->Write();
+  delete esdTree;
+  ef->Close();
+
+  delete vertexer;
+  delete rl;
+
+  return;
+}
+//--------------------------------------------------------------------------
+void BackToTPCInnerWall(AliESDEvent *event,AliESDEvent *eventTPC) {
+
+  Int_t ntracks = eventTPC->GetNumberOfTracks();
+  AliESDtrack *esdTrackTPC = 0;
+
+  // create relation between event and eventTPC
+  Int_t labelsTPC[100000000];
+  for(Int_t tr = 0; tr<ntracks; tr++) {
+    esdTrackTPC = (AliESDtrack*)event->GetTrack(tr);
+    labelsTPC[TMath::Abs(esdTrackTPC->GetLabel())] = tr;
+  }
+
+  ntracks = event->GetNumberOfTracks();
+  AliESDtrack *esdTrack = 0;
+  esdTrackTPC = 0;
+  Int_t indexTPC;
+
+  // loop on tracks
+  for(tr = 0; tr<ntracks; tr++) {
+    esdTrack = (AliESDtrack*)event->GetTrack(tr);
+    // set to kITSout the tracks that don't have kTPCin
+    // (they've been found by AliITStrackerSA)
+    if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) {
+      esdTrack->SetStatus(AliESDtrack::kITSout);
+      continue;
+    }
+
+    // skip tracks that don't have kITSin
+    if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
+
+    indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
+    esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
+
+    AliITStrackMI *itsTrack = 0;
+    try {
+      itsTrack = new AliITStrackMI(*esdTrackTPC);
+      esdTrack = 0;
+    }
+    catch (const Char_t *msg) {
+        Warning("ToTPCInnerWall",msg);
+        continue;
+    }
+    itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
+    esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
+
+    delete itsTrack;
+  } // end loop on tracks
+
+  return;
 }
+//--------------------------------------------------------------------------