[u/mrichter/AliRoot.git] / TPC / AliBarrelReconstruction.C

/****************************************************************************
 * This macro is supposed to do reconstruction in the barrel ALICE trackers *
 * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
 * It does the following steps (April 12, 2001):                            *
 *                   1) TPC cluster finding                                 *
 *                   2) TPC track finding                                   *
 *                   3) ITS cluster finding V2 (via fast points !)          *
 *                   4) ITS track finding V2                                *
 *                   5) ITS back track propagation V2                       *
 *                   6) TPC back track propagation                          *
 *                (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch)          *
 ****************************************************************************/

#ifndef __CINT__
  #include "alles.h"
  #include "iostream.h"
  #include "AliRun.h"
  #include "AliRunLoader.h"
  #include "AliLoader.h"
  #include "AliTPCLoader.h"
  #include "AliITSLoader.h"
  #include "AliMagF.h"
  #include "AliTPCtracker.h"
  #include "AliITS.h"
  #include "AliITSgeom.h"
  #include "AliITSRecPoint.h"
  #include "AliITSclusterV2.h"
  #include "AliITSsimulationFastPoints.h"
  #include "AliITStrackerV2.h"
#endif


Int_t TPCFindClusters( Int_t n);
Int_t TPCFindTracks(Int_t n);
Int_t TPCSortTracks(const Char_t *outname,  Int_t n);
Int_t TPCPropagateBack();

Int_t ITSFindClusters(Int_t n);
Int_t ITSFindTracks(const Char_t *inname2, Int_t n);
Int_t ITSPropagateBack();

const char* TPCtrkNameS= "TPC.TracksSorted.root";

class AliRunLoader;
class AliTPCLoader;
class AliTPCParam;
AliRunLoader *rl = 0x0;
AliTPCLoader *tpcl = 0x0;   
AliTPCParam  *param = 0x0;
Bool_t debug = kFALSE;
Int_t AliBarrelReconstruction(Int_t n=1) 
 {
 
   if (gAlice)
    {
      delete gAlice->GetRunLoader();
      delete gAlice;//if everything was OK here it is already NULL
      gAlice = 0x0;
    }
    
   rl = AliRunLoader::Open("galice.root");
   if (rl == 0x0)
    {
      cerr<<"Can not open session"<<endl;
      return 1;
    }
   
   if (rl->LoadgAlice())
    {
      cerr<<"Error occured while l"<<endl;
      return 1;
    }
   AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
   
   tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
   if (tpcl == 0x0)
    {
      cerr<<"Can not get TPC Loader"<<endl;
      return 1;
    }

// ********** Find TPC clusters *********** //
   if ( TPCFindClusters(n) ) 
    {
      cerr<<"Failed to get TPC clusters !\n";
      return 1;
    }
   
// ********** Find TPC tracks *********** //
    if (TPCFindTracks(n)) 
     {
      cerr<<"Failed to get TPC tracks !\n";
      return 1;
     }

//   cout<<"Stopping tracking on TPC\n";
// ********** Sort and label TPC tracks *********** //
   if (TPCSortTracks(TPCtrkNameS,n)) {
      cerr<<"Failed to sort TPC tracks !\n";
      return 1;
    } 

//   cout<<"Stopping tracking on TPC sorted T\n";

   
// ********** Find ITS clusters *********** //
   if (ITSFindClusters(n)) 
    {
      cerr<<"Failed to get ITS clusters !\n";
      return 1;
    }
   
// ********** Find ITS tracks *********** //

   if (ITSFindTracks(TPCtrkNameS,n)) 
    {
      cerr<<"Failed to get ITS tracks !\n";
      return 1;
    }
   //clsFile->Close();
   cout<<"Stopping on ITSFindTracks\n";
   delete rl;
   rl = 0x0;
   return 0;
   return 1;

// ********** Back propagation of the ITS tracks *********** //
   if ( ITSPropagateBack() ) {
      cerr<<"Failed to propagate back the ITS tracks !\n";
      return 1;
   }


// ********** Back propagation of the TPC tracks *********** //
   
   if (TPCPropagateBack()) {
      cerr<<"Failed to propagate back the TPC tracks !\n";
      return 1;
   }

   return 0;
}


Int_t TPCFindClusters(Int_t n) 
 {
   Int_t rc=0;
   const Char_t *name="TPCFindClusters";
   cerr<<'\n'<<name<<"...\n";
   rl->CdGAFile();
   param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
   if (!param) 
    {
     param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
     if (!param) 
      {
        cerr<<"TPC parameters have not been found !\n";
        return 1;
      }
    }
//   param->Dump();
   param->Update();

   gBenchmark->Start(name);

   AliTPCv2 tpc;
   tpc.SetParam(param);
   tpc.SetLoader(tpcl);
 
   tpcl->LoadDigits("read");
   tpcl->LoadRecPoints("recreate");
  

   //tpc.Digits2Clusters(out); //MI change
   for (Int_t i=0;i<n;i++)
    {
      printf("Processing event %d\n",i);
      tpc.Digits2Clusters(i);
     //     AliTPCclusterer::Digits2Clusters(dig, out, i);
    }
   
   tpcl->UnloadDigits();
   tpcl->UnloadRecPoints();
   gBenchmark->Stop(name);
   gBenchmark->Show(name);

   return rc;
}

Int_t TPCFindTracks(Int_t n) {

   Int_t rc=0;
   AliTPCtracker *tracker = 0x0;
   const Char_t *name="TPCFindTracks";
   rl->GetEvent(0);
   rl->CdGAFile();
   param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
   
   if (!param) 
    {
     param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
     if (!param) 
      {
        cerr<<"TPC parameters have not been found !\n";
        return 1;
      }
    }
//   param->Dump();
   param->Update();

   gBenchmark->Start(name);

   tpcl->LoadRecPoints("read");
   tpcl->LoadTracks("recreate");

   //AliTPCtracker *tracker=new AliTPCtracker(param);
   //rc=tracker->Clusters2Tracks(0,out);
   //delete tracker;
   
   for (Int_t i=0;i<n;i++){
     rl->GetEvent(i);
     printf("Processing event %d\n",i);
     
     tracker = new AliTPCtracker(param, i, (AliConfig::GetDefaultEventFolderName()).Data());
     //Int_t rc=
     tracker->Clusters2Tracks();
     delete tracker;
   }
   tpcl->UnloadRecPoints();
   tpcl->UnloadTracks();

   gBenchmark->Stop(name);
   gBenchmark->Show(name);

   return rc;
}


Int_t TPCSortTracks(const Char_t * outname,  Int_t eventn){
   Int_t rc=0;
   const Char_t *name="TPCSortTracks";

   cerr<<'\n'<<name<<"...\n";
   rl->CdGAFile();
   param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
   if (!param) 
    {
     param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
     if (!param) 
      {
        cerr<<"TPC parameters have not been found !\n";
        return 1;
      }
    }
   param->Update();

   gBenchmark->Start(name);
   
 
   AliRunLoader* rl2 = AliRunLoader::Open("galice.root","tmp");
   
   AliLoader* tpcl2 = (AliTPCLoader*)rl2->GetLoader("TPCLoader");
   cout<<"tpcl2->SetTracksFileName("<<outname<<");\n";
   tpcl2->SetTracksFileName(TString(outname));
   tpcl2->LoadTracks("recreate");
   
   // loop over events 
   for (Int_t event=0;event<eventn; event++)
    {
     
     rl->GetEvent(event);
     rl2->GetEvent(event);

     TObjArray tarray(10000);
     AliTPCtrack *iotrack=0;
     Int_t i;
     
     if (tpcl->TreeT() == 0x0)   tpcl->LoadTracks("read");
     TTree *tracktree=tpcl->TreeT();
     if (tracktree == 0x0)
      {
        cerr<<"Can not get TreeT for event "<<event<<endl;
        continue;
      }
     TBranch *tbranch=tracktree->GetBranch("tracks");
     Int_t nentr=(Int_t)tracktree->GetEntries();
     cout<<"IN Tracks nentr = "<<nentr<<endl;
     for (i=0; i<nentr; i++) 
      {
       iotrack=new AliTPCtrack;
       tbranch->SetAddress(&iotrack);
       tracktree->GetEvent(i);
       tarray.AddLast(iotrack);
      }   
     tarray.Sort();
     
     //assign thacks GEANT labels
     cout<<"Running Tracker\n";
     AliTPCtracker *tracker = new AliTPCtracker(param, event);

     cout<<"Load Sectors\n";

     tracker->LoadInnerSectors();
     tracker->LoadOuterSectors();

     cout<<"Cooking Labels\n";

     for (i=0; i<nentr; i++) 
      {
       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
       tracker->CookLabel(iotrack,0.1);
      }   
     cout<<"deleting tracker\n";
     delete tracker;
     //in->Close();
     //end of GEANT label assignment
     

     if (tpcl2->TreeT() == 0x0)   tpcl2->MakeTree("T");
     tracktree= tpcl2->TreeT();
     if (tracktree == 0x0)
      {
        cerr<<"Can not get TreeT for Sorted tracks\n";
        return 1;
      }
     tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
     for (i=0; i<nentr; i++) {
       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
       tracktree->Fill();
     }

     cout<<" Tracks File Name is "<<tpcl2->GetTracksFileName()<<endl;
     tpcl2->WriteTracks("OVERWRITE");
   }
  

   delete rl2;

   tpcl->UnloadTracks();
   tpcl->UnloadRecPoints(); 
   
   gBenchmark->Stop(name);
   gBenchmark->Show(name);

   return rc;
}

Int_t ITSFindClusters(Int_t n) 
 {
   Int_t rc=0;
   const Char_t *name="ITSFindClusters";
   cerr<<'\n'<<name<<"...\n";
   gBenchmark->Start(name);

   if (rl->GetEvent(0))
    {
      cerr<<"Problems\n";
      return 1;
    }
    
   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
   if (itsl == 0x0)
    {
      cerr<<"Can not get ITS Loader from Run Loader\n";
      return 1;
    } 
  
   gAlice=rl->GetAliRun();
   if (!gAlice) {
      cerr<<"Can't get gAlice !\n";
      return 1;
   }
   if (rl->TreeK() == 0x0) rl->LoadKinematics();
   if (rl->TreeE() == 0x0) rl->LoadHeader();
   
   itsl->LoadRawClusters("recreate");
   
   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
   if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
   AliITSgeom *geom=ITS->GetITSgeom();

   Int_t ev=0;
   for (ev = 0; ev<n; ev++)
   {

     rl->GetEvent(ev);
     TBranch *branch = 0x0;
     TTree *pTree = 0x0;
     
     pTree=itsl->TreeR();
     if (pTree == 0x0) 
      {
        itsl->LoadRecPoints("read");
        pTree=itsl->TreeR();
        if (pTree == 0x0) 
         {
           ::Error("AliBarrelReonstruction.C::ITSFindClusters",
                   "Can not get TreeR for event %d",ev);
           return 1;
         }
      }
     
     pTree=itsl->TreeR();
     branch=(pTree)?pTree->GetBranch("ITSRecPoints"):0x0;

     if (branch== 0x0) {
       //if not reconstructed ITS branch do reconstruction 
       if (debug) ::Info("AliBarrelReconstruction.C","Did not get ITSRecPoints from TreeR.");
       if (debug) ::Info("AliBarrelReconstruction.C","Making branch and Producing RecPoints");
       itsl->SetRecPointsFileOption("recreate");
       if (itsl->TreeR()==0x0) itsl->MakeTree("R");
       
       ITS->MakeBranch("RF");
       ITS->SetTreeAddress();
       itsl->LoadHits();

       //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
       for (Int_t i=0;i<3;i++) { 
         ITS->SetSimulationModel(i,new AliITSsimulationFastPoints()); 
       }
       Int_t nsignal=25;
       Int_t size=-1;
       Int_t bgr_ev=Int_t(ev/nsignal);
       
       ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
       //////////////////////////////////////////////////////////////////////////
       //MI comment  - in HitsToFast... they reset treeR to 0 
       //they overwrite full reconstructed event ???? ... so lets connect TreeR one more
       //time
       itsl->UnloadHits();
       ITS->SetTreeAddress();
     }

     TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);

     if (itsl->TreeC() == 0x0) itsl->MakeTree("C");
     TTree *cTree = itsl->TreeC();
     TBranch* b = cTree->GetBranch("Clusters");
     if (b == 0x0) 
       cTree->Branch("Clusters",&clusters);
     else b->SetAddress(&clusters);
     
      
     if (itsl->TreeR() == 0x0) itsl->LoadRecPoints();

     if (branch == 0x0)
      {//it means that we produced FastRecPoints above
       pTree=itsl->TreeR();
       if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
       branch=pTree->GetBranch("ITSRecPointsF");
       if (!branch) { cerr<<"Can't get Fast RecPoints branch named ITSRecPointsF  !\n"; return 1;}
      } 
     TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
     branch->SetAddress(&points);
     
     TClonesArray &cl=*clusters;
     Int_t nclusters=0;
     Int_t nentr=(Int_t)branch->GetEntries();
     AliITSgeom *geom=ITS->GetITSgeom();
     
     cout<<"\n\n Number of entries in TreeR = "<<nentr<<endl;
     for (Int_t i=0; i<nentr; i++) 
      {
       if (!branch->GetEvent(i)) 
        {
          cTree->Fill(); 
          continue;
        }
       
       Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
       Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
       Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
       Double_t yshift = x*rot[0] + y*rot[1];
       Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
       Int_t ncl=points->GetEntriesFast();
       nclusters+=ncl;

       if (debug) ::Info("AliBarrelReconstruction.C",
              "i=%d lay=%d lad=%d det=%d NRP=%d",
               i,   lay,   lad,   det,   ncl);

       for (Int_t j=0; j<ncl; j++) 
        { 
          Float_t lp[5];
          Int_t lab[4]; 

          AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
          lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
          lp[1]=p->GetZ()+zshift;
          lp[2]=p->GetSigmaX2();
          lp[3]=p->GetSigmaZ2();
          lp[4]=p->GetQ();
          lab[0]=p->GetLabel(0);
          lab[1]=p->GetLabel(1);
          lab[2]=p->GetLabel(2);
          lab[3]=ndet;
     
          Int_t label=lab[0];
          if (label<0) continue;
          TParticle *part=(TParticle*)rl->Stack()->Particle(label);
          label=-3;
          while (part->P() < 0.005) {
            Int_t m=part->GetFirstMother();
            if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
            label=m;
            part=(TParticle*)gAlice->Particle(label);
          }
          if (lab[1]<0) lab[1]=label;
          else if (lab[2]<0) lab[2]=label;
               else cerr<<"No empty labels !\n";
     
          new(cl[j]) AliITSclusterV2(lab,lp);
         }
       cTree->Fill(); clusters->Delete();
       points->Delete();
     }
     itsl->WriteRawClusters("OVERWRITE");
     
     cerr<<"Number of clusters: "<<nclusters<<endl;
     delete clusters; delete points;

   }
   itsl->UnloadRecPoints();
   itsl->UnloadRawClusters();
   rl->UnloadKinematics();

   gBenchmark->Stop(name);
   gBenchmark->Show(name);

   return rc;
}

Int_t ITSFindTracks(const Char_t *inname2, Int_t n) 
{

   Int_t rc=0;
   const Char_t *name="ITSFindTracks";
   cerr<<'\n'<<name<<"...\n";
   gBenchmark->Start(name);

   rl->GetEvent(0);

   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
   if ( (itsl == 0x0) || (tpcl == 0x0))
    {
      cerr<<"Can not get ITS Loader from Run Loader\n";
      return 1;
    } 
   
   tpcl->SetTracksFileName(TString(inname2));
   if(rl->GetAliRun() == 0x0) rl->LoadgAlice();
   
   AliITS* ITS = rl->GetAliRun()->GetModule("ITS");
   
   rl->CdGAFile();
   AliITSgeom *geom=ITS->GetITSgeom();
   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}

   itsl->LoadTracks("recreate");
   itsl->LoadRawClusters("read");
   for (Int_t i=0;i<n;i++)
    {
      rl->GetEvent(i);
      AliITStrackerV2 tracker(geom,i,AliConfig::GetDefaultEventFolderName());
      rc=tracker.Clusters2Tracks();
    }


   gBenchmark->Stop(name);
   gBenchmark->Show(name);

   return rc;
}

Int_t ITSPropagateBack() 
{
   
  Int_t rc=0;
  /*
   const Char_t *name="ITSPropagateBack";
   cerr<<'\n'<<name<<"...\n";
   gBenchmark->Start(name);

   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
   AliITStrackerV2 tracker(geom);

   TFile *out=TFile::Open(outname,"update");
   TFile *in =TFile::Open(inname);
   rc=tracker.PropagateBack(in,out);
   in->Close();
   out->Close();

   gBenchmark->Stop(name);
   gBenchmark->Show(name);
  */
   return rc;
}

Int_t TPCPropagateBack() 
{
   Int_t rc=0;
   const Char_t *name="TPCPropagateBack";
   cerr<<'\n'<<name<<"...\n";
   gBenchmark->Start(name);

   param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
   if (!param) 
    {
     param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
     if (!param) 
      {
        cerr<<"TPC parameters have not been found !\n";
        return 1;
      }
    }
//   param->Dump();
   param->Update();

   AliTPCtracker *tracker = new AliTPCtracker(param, AliConfig::GetDefaultEventFolderName());

//   TFile *out=TFile::Open(outname,"update");
//   TFile *in =TFile::Open(inname);

   rc=tracker->PropagateBack();
   delete tracker;

   gBenchmark->Stop(name);
   gBenchmark->Show(name);

   return rc;
}
Commit	Line	Data
b9de75e1	1	/****************************************************************************
	2	* This macro is supposed to do reconstruction in the barrel ALICE trackers *
	3	* (Make sure you have TPC digits and ITS hits before using this macro !!!) *
	4	* It does the following steps (April 12, 2001): *
	5	* 1) TPC cluster finding *
	6	* 2) TPC track finding *
	7	* 3) ITS cluster finding V2 (via fast points !) *
	8	* 4) ITS track finding V2 *
	9	* 5) ITS back track propagation V2 *
	10	* 6) TPC back track propagation *
	11	* (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch) *
	12	****************************************************************************/
	13
	14	#ifndef __CINT__
	15	#include "alles.h"
88cb7938	16	#include "iostream.h"
	17	#include "AliRun.h"
	18	#include "AliRunLoader.h"
	19	#include "AliLoader.h"
	20	#include "AliTPCLoader.h"
	21	#include "AliITSLoader.h"
b9de75e1	22	#include "AliMagF.h"
b9de75e1	23	#include "AliTPCtracker.h"
b9de75e1	24	#include "AliITS.h"
	25	#include "AliITSgeom.h"
	26	#include "AliITSRecPoint.h"
	27	#include "AliITSclusterV2.h"
	28	#include "AliITSsimulationFastPoints.h"
	29	#include "AliITStrackerV2.h"
	30	#endif
	31
b9de75e1	32
b9de75e1	33
88cb7938	34	Int_t TPCFindClusters( Int_t n);
	35	Int_t TPCFindTracks(Int_t n);
	36	Int_t TPCSortTracks(const Char_t *outname, Int_t n);
	37	Int_t TPCPropagateBack();
	38
	39	Int_t ITSFindClusters(Int_t n);
	40	Int_t ITSFindTracks(const Char_t *inname2, Int_t n);
	41	Int_t ITSPropagateBack();
	42
	43	const char* TPCtrkNameS= "TPC.TracksSorted.root";
	44
	45	class AliRunLoader;
	46	class AliTPCLoader;
	47	class AliTPCParam;
	48	AliRunLoader *rl = 0x0;
	49	AliTPCLoader *tpcl = 0x0;
	50	AliTPCParam *param = 0x0;
	51	Bool_t debug = kFALSE;
	52	Int_t AliBarrelReconstruction(Int_t n=1)
	53	{
	54
	55	if (gAlice)
	56	{
	57	delete gAlice->GetRunLoader();
	58	delete gAlice;//if everything was OK here it is already NULL
	59	gAlice = 0x0;
	60	}
743a19f2	61
88cb7938	62	rl = AliRunLoader::Open("galice.root");
	63	if (rl == 0x0)
	64	{
	65	cerr<<"Can not open session"<<endl;
	66	return 1;
	67	}
	68
	69	if (rl->LoadgAlice())
	70	{
	71	cerr<<"Error occured while l"<<endl;
	72	return 1;
	73	}
	74	AliKalmanTrack::SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
	75
	76	tpcl = (AliTPCLoader*)rl->GetLoader("TPCLoader");
	77	if (tpcl == 0x0)
	78	{
	79	cerr<<"Can not get TPC Loader"<<endl;
	80	return 1;
	81	}
	82
b9de75e1	83	// ******** Find TPC clusters ********* //
88cb7938	84	if ( TPCFindClusters(n) )
88cb7938	85	{
b9de75e1	86	cerr<<"Failed to get TPC clusters !\n";
b9de75e1	87	return 1;
88cb7938	88	}
88cb7938	89
b9de75e1	90	// ******** Find TPC tracks ********* //
88cb7938	91	if (TPCFindTracks(n))
88cb7938	92	{
b9de75e1	93	cerr<<"Failed to get TPC tracks !\n";
b9de75e1	94	return 1;
88cb7938	95	}
	96
	97	// cout<<"Stopping tracking on TPC\n";
b9de75e1	98	// ******** Sort and label TPC tracks ********* //
88cb7938	99	if (TPCSortTracks(TPCtrkNameS,n)) {
b9de75e1	100	cerr<<"Failed to sort TPC tracks !\n";
b9de75e1	101	return 1;
88cb7938	102	}
	103
	104	// cout<<"Stopping tracking on TPC sorted T\n";
	105
743a19f2	106
b9de75e1	107	// ******** Find ITS clusters ********* //
88cb7938	108	if (ITSFindClusters(n))
88cb7938	109	{
b9de75e1	110	cerr<<"Failed to get ITS clusters !\n";
b9de75e1	111	return 1;
88cb7938	112	}
743a19f2	113
b9de75e1	114	// ******** Find ITS tracks ********* //
88cb7938	115
	116	if (ITSFindTracks(TPCtrkNameS,n))
	117	{
b9de75e1	118	cerr<<"Failed to get ITS tracks !\n";
b9de75e1	119	return 1;
88cb7938	120	}
743a19f2	121	//clsFile->Close();
88cb7938	122	cout<<"Stopping on ITSFindTracks\n";
	123	delete rl;
	124	rl = 0x0;
	125	return 0;
743a19f2	126	return 1;
88cb7938	127
b9de75e1	128	// ******** Back propagation of the ITS tracks ********* //
88cb7938	129	if ( ITSPropagateBack() ) {
b9de75e1	130	cerr<<"Failed to propagate back the ITS tracks !\n";
	131	return 1;
	132	}
b9de75e1	133
	134
	135	// ******** Back propagation of the TPC tracks ********* //
88cb7938	136
88cb7938	137	if (TPCPropagateBack()) {
b9de75e1	138	cerr<<"Failed to propagate back the TPC tracks !\n";
	139	return 1;
	140	}
b9de75e1	141
	142	return 0;
	143	}
	144
	145
88cb7938	146	Int_t TPCFindClusters(Int_t n)
88cb7938	147	{
b9de75e1	148	Int_t rc=0;
	149	const Char_t *name="TPCFindClusters";
	150	cerr<<'\n'<<name<<"...\n";
88cb7938	151	rl->CdGAFile();
	152	param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
	153	if (!param)
	154	{
	155	param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
	156	if (!param)
	157	{
	158	cerr<<"TPC parameters have not been found !\n";
	159	return 1;
	160	}
	161	}
	162	// param->Dump();
	163	param->Update();
b9de75e1	164
88cb7938	165	gBenchmark->Start(name);
b9de75e1	166
b9de75e1	167	AliTPCv2 tpc;
b9de75e1	168	tpc.SetParam(param);
88cb7938	169	tpc.SetLoader(tpcl);
	170
	171	tpcl->LoadDigits("read");
	172	tpcl->LoadRecPoints("recreate");
	173
b9de75e1	174
743a19f2	175	//tpc.Digits2Clusters(out); //MI change
88cb7938	176	for (Int_t i=0;i<n;i++)
	177	{
	178	printf("Processing event %d\n",i);
	179	tpc.Digits2Clusters(i);
	180	// AliTPCclusterer::Digits2Clusters(dig, out, i);
	181	}
	182
	183	tpcl->UnloadDigits();
	184	tpcl->UnloadRecPoints();
b9de75e1	185	gBenchmark->Stop(name);
	186	gBenchmark->Show(name);
	187
	188	return rc;
	189	}
	190
88cb7938	191	Int_t TPCFindTracks(Int_t n) {
88cb7938	192
b9de75e1	193	Int_t rc=0;
88cb7938	194	AliTPCtracker *tracker = 0x0;
b9de75e1	195	const Char_t *name="TPCFindTracks";
88cb7938	196	rl->GetEvent(0);
	197	rl->CdGAFile();
	198	param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
	199
	200	if (!param)
	201	{
	202	param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
	203	if (!param)
	204	{
	205	cerr<<"TPC parameters have not been found !\n";
	206	return 1;
	207	}
	208	}
	209	// param->Dump();
	210	param->Update();
	211
b9de75e1	212	gBenchmark->Start(name);
88cb7938	213
	214	tpcl->LoadRecPoints("read");
	215	tpcl->LoadTracks("recreate");
b9de75e1	216
743a19f2	217	//AliTPCtracker *tracker=new AliTPCtracker(param);
	218	//rc=tracker->Clusters2Tracks(0,out);
	219	//delete tracker;
88cb7938	220
743a19f2	221	for (Int_t i=0;i<n;i++){
88cb7938	222	rl->GetEvent(i);
743a19f2	223	printf("Processing event %d\n",i);
88cb7938	224
f5a857b2	225	tracker = new AliTPCtracker(param, i, (AliConfig::GetDefaultEventFolderName()).Data());
743a19f2	226	//Int_t rc=
88cb7938	227	tracker->Clusters2Tracks();
743a19f2	228	delete tracker;
743a19f2	229	}
88cb7938	230	tpcl->UnloadRecPoints();
88cb7938	231	tpcl->UnloadTracks();
743a19f2	232
b9de75e1	233	gBenchmark->Stop(name);
	234	gBenchmark->Show(name);
	235
	236	return rc;
	237	}
743a19f2	238
743a19f2	239
88cb7938	240	Int_t TPCSortTracks(const Char_t * outname, Int_t eventn){
	241	Int_t rc=0;
	242	const Char_t *name="TPCSortTracks";
b9de75e1	243
88cb7938	244	cerr<<'\n'<<name<<"...\n";
	245	rl->CdGAFile();
	246	param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
	247	if (!param)
	248	{
	249	param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
	250	if (!param)
	251	{
	252	cerr<<"TPC parameters have not been found !\n";
	253	return 1;
	254	}
	255	}
	256	param->Update();
743a19f2	257
88cb7938	258	gBenchmark->Start(name);
	259
	260
	261	AliRunLoader* rl2 = AliRunLoader::Open("galice.root","tmp");
	262
	263	AliLoader* tpcl2 = (AliTPCLoader*)rl2->GetLoader("TPCLoader");
	264	cout<<"tpcl2->SetTracksFileName("<<outname<<");\n";
	265	tpcl2->SetTracksFileName(TString(outname));
	266	tpcl2->LoadTracks("recreate");
	267
743a19f2	268	// loop over events
88cb7938	269	for (Int_t event=0;event<eventn; event++)
	270	{
	271
	272	rl->GetEvent(event);
	273	rl2->GetEvent(event);
	274
743a19f2	275	TObjArray tarray(10000);
	276	AliTPCtrack *iotrack=0;
	277	Int_t i;
88cb7938	278
	279	if (tpcl->TreeT() == 0x0) tpcl->LoadTracks("read");
	280	TTree *tracktree=tpcl->TreeT();
	281	if (tracktree == 0x0)
	282	{
	283	cerr<<"Can not get TreeT for event "<<event<<endl;
	284	continue;
	285	}
743a19f2	286	TBranch *tbranch=tracktree->GetBranch("tracks");
743a19f2	287	Int_t nentr=(Int_t)tracktree->GetEntries();
88cb7938	288	cout<<"IN Tracks nentr = "<<nentr<<endl;
	289	for (i=0; i<nentr; i++)
	290	{
b9de75e1	291	iotrack=new AliTPCtrack;
	292	tbranch->SetAddress(&iotrack);
	293	tracktree->GetEvent(i);
	294	tarray.AddLast(iotrack);
88cb7938	295	}
743a19f2	296	tarray.Sort();
743a19f2	297
743a19f2	298	//assign thacks GEANT labels
88cb7938	299	cout<<"Running Tracker\n";
	300	AliTPCtracker *tracker = new AliTPCtracker(param, event);
	301
	302	cout<<"Load Sectors\n";
	303
743a19f2	304	tracker->LoadInnerSectors();
743a19f2	305	tracker->LoadOuterSectors();
88cb7938	306
	307	cout<<"Cooking Labels\n";
	308
	309	for (i=0; i<nentr; i++)
	310	{
b9de75e1	311	iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
b9de75e1	312	tracker->CookLabel(iotrack,0.1);
88cb7938	313	}
88cb7938	314	cout<<"deleting tracker\n";
743a19f2	315	delete tracker;
	316	//in->Close();
	317	//end of GEANT label assignment
	318
88cb7938	319
	320	if (tpcl2->TreeT() == 0x0) tpcl2->MakeTree("T");
	321	tracktree= tpcl2->TreeT();
	322	if (tracktree == 0x0)
	323	{
	324	cerr<<"Can not get TreeT for Sorted tracks\n";
	325	return 1;
	326	}
743a19f2	327	tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
743a19f2	328	for (i=0; i<nentr; i++) {
b9de75e1	329	iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
b9de75e1	330	tracktree->Fill();
743a19f2	331	}
88cb7938	332
	333	cout<<" Tracks File Name is "<<tpcl2->GetTracksFileName()<<endl;
	334	tpcl2->WriteTracks("OVERWRITE");
b9de75e1	335	}
88cb7938	336
743a19f2	337
88cb7938	338	delete rl2;
b9de75e1	339
88cb7938	340	tpcl->UnloadTracks();
	341	tpcl->UnloadRecPoints();
	342
b9de75e1	343	gBenchmark->Stop(name);
	344	gBenchmark->Show(name);
	345
	346	return rc;
	347	}
	348
88cb7938	349	Int_t ITSFindClusters(Int_t n)
88cb7938	350	{
b9de75e1	351	Int_t rc=0;
	352	const Char_t *name="ITSFindClusters";
	353	cerr<<'\n'<<name<<"...\n";
	354	gBenchmark->Start(name);
b9de75e1	355
88cb7938	356	if (rl->GetEvent(0))
	357	{
	358	cerr<<"Problems\n";
	359	return 1;
	360	}
	361
	362	AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
	363	if (itsl == 0x0)
	364	{
	365	cerr<<"Can not get ITS Loader from Run Loader\n";
	366	return 1;
	367	}
	368
	369	gAlice=rl->GetAliRun();
	370	if (!gAlice) {
b9de75e1	371	cerr<<"Can't get gAlice !\n";
	372	return 1;
	373	}
88cb7938	374	if (rl->TreeK() == 0x0) rl->LoadKinematics();
	375	if (rl->TreeE() == 0x0) rl->LoadHeader();
	376
	377	itsl->LoadRawClusters("recreate");
	378
b9de75e1	379	AliITS ITS = (AliITS)gAlice->GetModule("ITS");
b9de75e1	380	if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
b9de75e1	381	AliITSgeom *geom=ITS->GetITSgeom();
88cb7938	382
743a19f2	383	Int_t ev=0;
88cb7938	384	for (ev = 0; ev<n; ev++)
88cb7938	385	{
743a19f2	386
88cb7938	387	rl->GetEvent(ev);
	388	TBranch *branch = 0x0;
	389	TTree *pTree = 0x0;
743a19f2	390
88cb7938	391	pTree=itsl->TreeR();
	392	if (pTree == 0x0)
	393	{
	394	itsl->LoadRecPoints("read");
	395	pTree=itsl->TreeR();
	396	if (pTree == 0x0)
	397	{
	398	::Error("AliBarrelReonstruction.C::ITSFindClusters",
	399	"Can not get TreeR for event %d",ev);
	400	return 1;
	401	}
	402	}
	403
	404	pTree=itsl->TreeR();
	405	branch=(pTree)?pTree->GetBranch("ITSRecPoints"):0x0;
	406
	407	if (branch== 0x0) {
743a19f2	408	//if not reconstructed ITS branch do reconstruction
88cb7938	409	if (debug) ::Info("AliBarrelReconstruction.C","Did not get ITSRecPoints from TreeR.");
	410	if (debug) ::Info("AliBarrelReconstruction.C","Making branch and Producing RecPoints");
	411	itsl->SetRecPointsFileOption("recreate");
	412	if (itsl->TreeR()==0x0) itsl->MakeTree("R");
	413
	414	ITS->MakeBranch("RF");
	415	ITS->SetTreeAddress();
	416	itsl->LoadHits();
	417
743a19f2	418	//////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
8d9c76b5	419	for (Int_t i=0;i<3;i++) {
	420	ITS->SetSimulationModel(i,new AliITSsimulationFastPoints());
	421	}
743a19f2	422	Int_t nsignal=25;
	423	Int_t size=-1;
	424	Int_t bgr_ev=Int_t(ev/nsignal);
88cb7938	425
743a19f2	426	ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
743a19f2	427	//////////////////////////////////////////////////////////////////////////
88cb7938	428	//MI comment - in HitsToFast... they reset treeR to 0
743a19f2	429	//they overwrite full reconstructed event ???? ... so lets connect TreeR one more
743a19f2	430	//time
88cb7938	431	itsl->UnloadHits();
88cb7938	432	ITS->SetTreeAddress();
743a19f2	433	}
743a19f2	434
743a19f2	435	TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
88cb7938	436
	437	if (itsl->TreeC() == 0x0) itsl->MakeTree("C");
	438	TTree *cTree = itsl->TreeC();
	439	TBranch* b = cTree->GetBranch("Clusters");
	440	if (b == 0x0)
	441	cTree->Branch("Clusters",&clusters);
	442	else b->SetAddress(&clusters);
743a19f2	443
88cb7938	444
	445	if (itsl->TreeR() == 0x0) itsl->LoadRecPoints();
	446
	447	if (branch == 0x0)
	448	{//it means that we produced FastRecPoints above
	449	pTree=itsl->TreeR();
	450	if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
	451	branch=pTree->GetBranch("ITSRecPointsF");
	452	if (!branch) { cerr<<"Can't get Fast RecPoints branch named ITSRecPointsF !\n"; return 1;}
	453	}
743a19f2	454	TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
	455	branch->SetAddress(&points);
	456
	457	TClonesArray &cl=*clusters;
	458	Int_t nclusters=0;
88cb7938	459	Int_t nentr=(Int_t)branch->GetEntries();
743a19f2	460	AliITSgeom *geom=ITS->GetITSgeom();
88cb7938	461
	462	cout<<"\n\n Number of entries in TreeR = "<<nentr<<endl;
	463	for (Int_t i=0; i<nentr; i++)
	464	{
	465	if (!branch->GetEvent(i))
	466	{
	467	cTree->Fill();
	468	continue;
	469	}
	470
b9de75e1	471	Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
	472	Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift);
	473	Double_t rot[9]; geom->GetRotMatrix(lay,lad,det,rot);
	474	Double_t yshift = xrot[0] + yrot[1];
	475	Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
	476	Int_t ncl=points->GetEntriesFast();
	477	nclusters+=ncl;
88cb7938	478
	479	if (debug) ::Info("AliBarrelReconstruction.C",
	480	"i=%d lay=%d lad=%d det=%d NRP=%d",
	481	i, lay, lad, det, ncl);
	482
	483	for (Int_t j=0; j<ncl; j++)
	484	{
	485	Float_t lp[5];
	486	Int_t lab[4];
	487
	488	AliITSRecPoint p=(AliITSRecPoint)points->UncheckedAt(j);
	489	lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
	490	lp[1]=p->GetZ()+zshift;
	491	lp[2]=p->GetSigmaX2();
	492	lp[3]=p->GetSigmaZ2();
	493	lp[4]=p->GetQ();
	494	lab[0]=p->GetLabel(0);
	495	lab[1]=p->GetLabel(1);
	496	lab[2]=p->GetLabel(2);
	497	lab[3]=ndet;
	498
	499	Int_t label=lab[0];
	500	if (label<0) continue;
	501	TParticle part=(TParticle)rl->Stack()->Particle(label);
	502	label=-3;
	503	while (part->P() < 0.005) {
	504	Int_t m=part->GetFirstMother();
	505	if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
	506	label=m;
	507	part=(TParticle*)gAlice->Particle(label);
	508	}
	509	if (lab[1]<0) lab[1]=label;
	510	else if (lab[2]<0) lab[2]=label;
	511	else cerr<<"No empty labels !\n";
	512
	513	new(cl[j]) AliITSclusterV2(lab,lp);
	514	}
b9de75e1	515	cTree->Fill(); clusters->Delete();
b9de75e1	516	points->Delete();
743a19f2	517	}
88cb7938	518	itsl->WriteRawClusters("OVERWRITE");
88cb7938	519
743a19f2	520	cerr<<"Number of clusters: "<<nclusters<<endl;
88cb7938	521	delete clusters; delete points;
b9de75e1	522
743a19f2	523	}
88cb7938	524	itsl->UnloadRecPoints();
	525	itsl->UnloadRawClusters();
	526	rl->UnloadKinematics();
b9de75e1	527
b9de75e1	528	gBenchmark->Stop(name);
	529	gBenchmark->Show(name);
	530
	531	return rc;
	532	}
	533
88cb7938	534	Int_t ITSFindTracks(const Char_t *inname2, Int_t n)
	535	{
	536
b9de75e1	537	Int_t rc=0;
	538	const Char_t *name="ITSFindTracks";
	539	cerr<<'\n'<<name<<"...\n";
	540	gBenchmark->Start(name);
	541
88cb7938	542	rl->GetEvent(0);
b9de75e1	543
88cb7938	544	AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
	545	if ( (itsl == 0x0) \|\| (tpcl == 0x0))
	546	{
	547	cerr<<"Can not get ITS Loader from Run Loader\n";
	548	return 1;
	549	}
	550
	551	tpcl->SetTracksFileName(TString(inname2));
	552	if(rl->GetAliRun() == 0x0) rl->LoadgAlice();
	553
	554	AliITS* ITS = rl->GetAliRun()->GetModule("ITS");
	555
	556	rl->CdGAFile();
	557	AliITSgeom *geom=ITS->GetITSgeom();
743a19f2	558	if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
743a19f2	559
88cb7938	560	itsl->LoadTracks("recreate");
	561	itsl->LoadRawClusters("read");
	562	for (Int_t i=0;i<n;i++)
	563	{
	564	rl->GetEvent(i);
f5a857b2	565	AliITStrackerV2 tracker(geom,i,AliConfig::GetDefaultEventFolderName());
88cb7938	566	rc=tracker.Clusters2Tracks();
88cb7938	567	}
743a19f2	568
b9de75e1	569
	570	gBenchmark->Stop(name);
	571	gBenchmark->Show(name);
	572
	573	return rc;
	574	}
	575
88cb7938	576	Int_t ITSPropagateBack()
88cb7938	577	{
743a19f2	578
	579	Int_t rc=0;
	580	/*
b9de75e1	581	const Char_t *name="ITSPropagateBack";
	582	cerr<<'\n'<<name<<"...\n";
	583	gBenchmark->Start(name);
	584
	585	AliITSgeom geom=(AliITSgeom)gFile->Get("AliITSgeom");
	586	if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
	587	AliITStrackerV2 tracker(geom);
	588
	589	TFile *out=TFile::Open(outname,"update");
	590	TFile *in =TFile::Open(inname);
	591	rc=tracker.PropagateBack(in,out);
	592	in->Close();
	593	out->Close();
	594
	595	gBenchmark->Stop(name);
	596	gBenchmark->Show(name);
743a19f2	597	*/
b9de75e1	598	return rc;
	599	}
	600
88cb7938	601	Int_t TPCPropagateBack()
88cb7938	602	{
b9de75e1	603	Int_t rc=0;
	604	const Char_t *name="TPCPropagateBack";
	605	cerr<<'\n'<<name<<"...\n";
	606	gBenchmark->Start(name);
	607
88cb7938	608	param=(AliTPCParam *)gDirectory->Get("75x40_100x60");
	609	if (!param)
	610	{
	611	param=(AliTPCParam *)gDirectory->Get("75x40_100x60_150x60");
	612	if (!param)
	613	{
	614	cerr<<"TPC parameters have not been found !\n";
	615	return 1;
	616	}
	617	}
	618	// param->Dump();
	619	param->Update();
b9de75e1	620
f5a857b2	621	AliTPCtracker *tracker = new AliTPCtracker(param, AliConfig::GetDefaultEventFolderName());
88cb7938	622
	623	// TFile *out=TFile::Open(outname,"update");
	624	// TFile *in =TFile::Open(inname);
	625
	626	rc=tracker->PropagateBack();
b9de75e1	627	delete tracker;
b9de75e1	628
	629	gBenchmark->Stop(name);
	630	gBenchmark->Show(name);
	631
	632	return rc;
	633	}