]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliESDtest.C
Reconstruction and PID using transition radiation photons: first implementation ...
[u/mrichter/AliRoot.git] / STEER / AliESDtest.C
index da6ddf1dc6969590631d43f7074ea87dce105b31..ad9c44ff2024629859fcb951d543168f214fcccf 100644 (file)
   #include "TFile.h"
   #include "TSystem.h"
   #include "TStopwatch.h"
-  #include "TGeant3.h"
+  #include "TArrayF.h"
 
   #include "AliMagF.h"
   #include "AliRun.h"
   #include "AliRunLoader.h"
   #include "AliLoader.h"
+  #include "AliHeader.h"
+  #include "AliGenEventHeader.h"
 
   #include "AliESD.h"
   #include "AliESDpid.h"
   #include "AliTRDPartID.h"
 
   #include "AliTOFpidESD.h"
+  #include "AliTOF.h"
+  #include "AliTOFGeometry.h"
 #endif
 
-extern TSystem *gSystem;
 extern AliRun *gAlice;
 extern TFile *gFile;
 
-Int_t AliESDtest(Int_t nev=1) {
+Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
 
 /**** Initialization of the NewIO *******/
 
@@ -63,9 +66,6 @@ Int_t AliESDtest(Int_t nev=1) {
       gAlice=0;
    }
 
-   gSystem->Load("libgeant321");     // needed for the PID in TOF 
-   new TGeant3("");                  // must be re-done !
-
    AliRunLoader *rl = AliRunLoader::Open("galice.root");
    if (rl == 0x0) {
       cerr<<"Can not open session"<<endl;
@@ -109,11 +109,9 @@ Int_t AliESDtest(Int_t nev=1) {
 
    //An instance of the ITS tracker
    AliITStrackerV2 itsTracker(geom);
-   {Double_t xyz[]={0.,0.,0.}, ers[]={0.005, 0.005, 0.010};
-   itsTracker.SetVertex(xyz,ers);}
    
    //An instance of the ITS PID maker
-   Double_t parITS[]={34.,0.15,10.};
+   Double_t parITS[]={35.5,0.11,10.};
    AliITSpidESD itsPID(parITS);
 
    //An instance of the V0 finder
@@ -158,7 +156,7 @@ Int_t AliESDtest(Int_t nev=1) {
    AliTPCtracker tpcTracker(par);
 
    //An instance of the TPC PID maker
-   Double_t parTPC[]={47.,0.10,10.};
+   Double_t parTPC[]={45.0,0.08,10.};
    AliTPCpidESD tpcPID(parTPC);
 
 
@@ -191,6 +189,16 @@ Int_t AliESDtest(Int_t nev=1) {
 
 
 /**** The TOF corner ********************/
+   AliTOF *dTOF = (AliTOF*)gAlice->GetDetector("TOF");
+   if (!dTOF) {
+      cerr<<"AliESDtest.C : Can not find the TOF detector !"<<endl;
+      return 4;
+   }
+   AliTOFGeometry *tofGeo = dTOF->GetGeometry();
+   if (!tofGeo) {
+      cerr<<"AliESDtest.C : Can not find the TOF geometry !"<<endl;
+      return 4;
+   }
 
    AliLoader* tofl = rl->GetLoader("TOFLoader");
    if (tofl == 0x0) {
@@ -199,30 +207,46 @@ Int_t AliESDtest(Int_t nev=1) {
    }
    tofl->LoadDigits("read");
 
-
    //Instance of the TOF PID maker
    Double_t parTOF[]={130.,5.};
-   AliTOFpidESD tofPID(parTOF);
+   AliTOFtracker tofTracker(tofGeo,parTOF) ;
 
 
    //rl->UnloadgAlice();
 
 
+   TFile *bf=TFile::Open("AliESDcheck.root","RECREATE");
+   if (!bf || !bf->IsOpen()) {
+      cerr<<"Can't open AliESDcheck.root !\n"; return 1;
+   }
    TFile *ef=TFile::Open("AliESDs.root","RECREATE");
-   if (!ef->IsOpen()) {cerr<<"Can't AliESDs.root !\n"; return 1;}
+   if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return 2;}
 
    TStopwatch timer;
    Int_t rc=0;
    if (nev>rl->GetNumberOfEvents()) nev=rl->GetNumberOfEvents();
    //The loop over events
    for (Int_t i=0; i<nev; i++) {
+     Char_t ename[100]; 
+
      cerr<<"\n\nProcessing event number : "<<i<<endl;
      AliESD *event=new AliESD(); 
-
+     event->SetRunNumber(run);
+     event->SetEventNumber(i);
+     event->SetMagneticField(gAlice->Field()->SolenoidField());
+     
      rl->GetEvent(i);
  
+//***** Primary vertex reconstruction (MC vertex position, for the moment)
+     TArrayF v(3);     
+     rl->GetHeader()->GenEventHeader()->PrimaryVertex(v);
+     Double_t vtx[3]={v[0],v[1],v[2]};
+     Double_t cvtx[3]={0.005,0.005,0.010};
+     AliESDVertex vertex(vtx,cvtx);
+     event->SetVertex(&vertex);
 
 //***** Initial path towards the primary vertex
+     tpcTracker.SetVertex(vtx,cvtx);
      TTree *tpcTree=tpcl->TreeR();
      if (!tpcTree) {
         cerr<<"Can't get the TPC cluster tree !\n";
@@ -230,7 +254,10 @@ Int_t AliESDtest(Int_t nev=1) {
      }     
      tpcTracker.LoadClusters(tpcTree);
      rc+=tpcTracker.Clusters2Tracks(event);
+     tpcPID.MakePID(event);                 // preliminary PID
+     AliESDpid::MakePID(event);             // for the ITS tracker
 
+     itsTracker.SetVertex(vtx,cvtx);
      TTree *itsTree=itsl->TreeR();
      if (!itsTree) {
         cerr<<"Can't get the ITS cluster tree !\n";
@@ -239,13 +266,16 @@ Int_t AliESDtest(Int_t nev=1) {
      itsTracker.LoadClusters(itsTree);
      rc+=itsTracker.Clusters2Tracks(event);
 
+       //checkpoint
+       bf->cd();
+       sprintf(ename,"in%d",i);
+       event->Write(ename); bf->Flush();
+       ef->cd();
 
 //***** Back propagation towards the outer barrel detectors
      rc+=itsTracker.PropagateBack(event); 
-     //itsPID.MakePID(event);
      
      rc+=tpcTracker.PropagateBack(event);
-     tpcPID.MakePID(event);
 
      TTree *trdTree=trdl->TreeR();
      if (!trdTree) {
@@ -254,27 +284,30 @@ Int_t AliESDtest(Int_t nev=1) {
      } 
      trdTracker.LoadClusters(trdTree);
      rc+=trdTracker.PropagateBack(event);
+
 /*
      for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) {
        AliESDtrack* track = event->GetTrack(iTrack);
        trdPID->MakePID(track);
      }
 */
+
      TTree *tofTree=tofl->TreeD();
      if (!tofTree) {
         cerr<<"Can't get the TOF cluster tree !\n";
         return 4;
      } 
-     tofPID.LoadClusters(tofTree);
-     tofPID.MakePID(event);
-     tofPID.UnloadClusters();
-
-
-
-//***** Here is the combined PID
-     AliESDpid::MakePID(event);
-
-
+     tofTracker.LoadClusters(tofTree);
+     rc+=tofTracker.PropagateBack(event);
+     tofTracker.UnloadClusters();
+
+       //checkpoint
+     bf->cd();
+     strcat(ename,";*"); bf->Delete(ename);
+     sprintf(ename,"out%d",i);
+     event->Write(ename); bf->Flush();
+     ef->cd();
+     
 
 //***** Now the final refit at the primary vertex...
      rc+=trdTracker.RefitInward(event);
@@ -282,27 +315,42 @@ Int_t AliESDtest(Int_t nev=1) {
 
      rc+=tpcTracker.RefitInward(event);
      tpcTracker.UnloadClusters();
+     tpcPID.MakePID(event);
 
      rc+=itsTracker.RefitInward(event); 
      itsTracker.UnloadClusters();
+     itsPID.MakePID(event);
+
 
+//***** Here is the combined PID
+     AliESDpid::MakePID(event);
 
 
+       //checkpoint
+     bf->cd();
+     strcat(ename,";*"); bf->Delete(ename);
+     sprintf(ename,"refit%d",i);
+     event->Write(ename); bf->Flush();
+     ef->cd();
+     
+     bf->Close();
+     
 //***** Hyperon reconstruction 
+     vtxer.SetVertex(vtx);
      rc+=vtxer.Tracks2V0vertices(event);            // V0 finding
      rc+=cvtxer.V0sTracks2CascadeVertices(event);   // cascade finding
 
 
-
 //***** Some final manipulations with this event 
      if (rc==0) {
-        Char_t ename[100]; 
         sprintf(ename,"%d",i);
        ef->cd();
-        if (!event->Write(ename)) rc++;
+        if (!event->Write(ename)) rc++; ef->Flush();
+        bf=TFile::Open("AliESDcheck.root","RECREATE");
      } 
      if (rc) {
         cerr<<"Something bad happened...\n";
+        bf=TFile::Open("AliESDcheck.root","UPDATE");
      }
      delete event;
    }
@@ -313,6 +361,7 @@ Int_t AliESDtest(Int_t nev=1) {
    delete par;
 
    ef->Close();
+   bf->Close();
 
    delete rl;