Reconstruction and PID using transition radiation photons: first implementation ...

[u/mrichter/AliRoot.git] / STEER / AliESDtest.C

diff --git a/STEER/AliESDtest.C b/STEER/AliESDtest.C
index d4c60c3926114253a8b49349b51492f6a32fb066..ad9c44ff2024629859fcb951d543168f214fcccf 100644 (file)
--- a/STEER/AliESDtest.C
+++ b/STEER/AliESDtest.C
@@ -20,7 +20,6 @@
   #include "TFile.h"
   #include "TSystem.h"
   #include "TStopwatch.h"
-  #include "TGeant3.h"
   #include "TArrayF.h"
 
   #include "AliMagF.h"
@@ -50,9 +49,10 @@
   #include "AliTRDPartID.h"
 
   #include "AliTOFpidESD.h"
+  #include "AliTOF.h"
+  #include "AliTOFGeometry.h"
 #endif
 
-extern TSystem *gSystem;
 extern AliRun *gAlice;
 extern TFile *gFile;
 
@@ -66,9 +66,6 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
       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;
@@ -114,7 +111,7 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
    AliITStrackerV2 itsTracker(geom);
    
    //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
@@ -159,7 +156,7 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
    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);
 
 
@@ -192,6 +189,16 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
 
 
 /**** 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) {
@@ -200,41 +207,43 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
    }
    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[6]={
-       0.005,
-       0.000, 0.005,
-       0.000, 0.000, 0.010
-     };
-     event->SetVertex(vtx,cvtx);
-     cvtx[1]=cvtx[0]; cvtx[2]=cvtx[5]; //trackers use only the diag.elements
+     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);
@@ -245,6 +254,8 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
      }     
      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();
@@ -255,13 +266,16 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
      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) {
@@ -270,27 +284,30 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
      } 
      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);
@@ -298,28 +315,42 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
 
      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;
    }
@@ -330,6 +361,7 @@ Int_t AliESDtest(Int_t nev=1,Int_t run=0) {
    delete par;
 
    ef->Close();
+   bf->Close();
 
    delete rl;