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

void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
  //
  // Macro to create AliESDs.root using parametrized TPC tracking
  // and AliITStrackerSA (MI + SA)
  //
  // Input files:
  // - galice.root
  // - Kinematics.root
  // - TrackReferences.root
  // - ITS.RecPoints.root (use AliRecontruction class)
  // - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
  //
  // A. Dainese - INFN Legnaro
  //

  Int_t  collcode = 1; // pp collisions
  Bool_t useMeanVtx = kFALSE;
  
  TGeoManager::Import("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");
  
  AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
  if (!dITS) {
    cerr<<"Cannot find the ITS detector !"<<endl;
    return;
  }
  AliITSgeom *geom = dITS->GetITSgeom();
  
  // Instance of the ITS tracker
  AliITStrackerSA itsTracker(geom);
  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();
  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");
  AliESD *event=0; AliESD *eventTPCin=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 *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;
    AliESD *event = new AliESD(); 
    event->SetRunNumber(gAlice->GetRunNumber());
    event->SetEventNumber(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(); delete event;
      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(); delete event;
      continue;
    }

    // make a copy of the ESD at this stage
    eventTPCin = event;

    //***** ITS tracking
    itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
    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;
    }

    // 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();
    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;}

  //Write the tree and close everything
  esdTree->Write();
  delete esdTree;
  ef->Close();

  delete vertexer;
  delete rl;

  return;
}
//--------------------------------------------------------------------------
void BackToTPCInnerWall(AliESD *event,AliESD *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;
}
//--------------------------------------------------------------------------
Commit	Line	Data
a654e257	1	void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
056891c0	2	//
a654e257	3	// Macro to create AliESDs.root using parametrized TPC tracking
970d08c3	4	// and AliITStrackerSA (MI + SA)
056891c0	5	//
a654e257	6	// Input files:
	7	// - galice.root
	8	// - Kinematics.root
	9	// - TrackReferences.root
	10	// - ITS.RecPoints.root (use AliRecontruction class)
	11	// - ITS.Vertex.root (use $ALICE_ROOT/ITS/AliITSVertexerZTest.C)
e0cf7c8a	12	//
970d08c3	13	// A. Dainese - INFN Legnaro
e0cf7c8a	14	//
056891c0	15
970d08c3	16	Int_t collcode = 1; // pp collisions
225a9242	17	Bool_t useMeanVtx = kFALSE;
970d08c3	18
2fdadb1c	19	TGeoManager::Import("geometry.root");
970d08c3	20
	21	if (gAlice) {
	22	delete gAlice->GetRunLoader();
	23	delete gAlice;
	24	gAlice=0;
	25	}
	26	AliRunLoader *rl = AliRunLoader::Open("galice.root");
	27	if (rl == 0x0) {
	28	cerr<<"Can not open session"<<endl;
	29	return;
	30	}
	31	Int_t retval = rl->LoadgAlice();
	32	if (retval) {
	33	cerr<<"LoadgAlice returned error"<<endl;
	34	delete rl;
	35	return;
	36	}
	37	retval = rl->LoadHeader();
	38	if (retval) {
	39	cerr<<"LoadHeader returned error"<<endl;
	40	delete rl;
	41	return;
	42	}
	43	gAlice=rl->GetAliRun();
	44
	45
	46
	47	// Get field from galice.root
	48	AliMagF fiel = (AliMagF)gAlice->Field();
	49	Double_t fieval=TMath::Abs((Double_t)fiel->SolenoidField()/10.);
	50	// Set the conversion constant between curvature and Pt
	51	AliTracker::SetFieldMap(fiel,kTRUE);
	52
	53	/** The TPC corner ******************/
	54
	55	AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
	56	tpcTrackerPar.Init();
	57
	58	//**** Switch on the PID class (mandatory!)
	59	AliPID pid;
	60
	61	/** The ITS corner ******************/
	62
	63	AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
	64	if (itsl == 0x0) {
	65	cerr<<"Cannot get the ITS loader"<<endl;
	66	return;
	67	}
	68	itsl->LoadRecPoints("read");
	69
	70	AliITS dITS = (AliITS)gAlice->GetDetector("ITS");
	71	if (!dITS) {
	72	cerr<<"Cannot find the ITS detector !"<<endl;
	73	return;
	74	}
	75	AliITSgeom *geom = dITS->GetITSgeom();
	76
	77	// Instance of the ITS tracker
	78	AliITStrackerSA itsTracker(geom);
	79	Int_t ITSclusters[6] = {1,1,1,1,1,1};
	80	itsTracker.SetLayersNotToSkip(ITSclusters);
	81
	82
	83	// Primary vertex reconstruction in pp
84	AliESDVertex *initVertex = 0;
85	TFile *invtx = new TFile("AliESDVertexMean.root");
86	if(collcode==1 && useMeanVtx) {
87	// open the mean vertex
88	initVertex = (AliESDVertex*)invtx->Get("vtxmean");
89	} else {
90	Double_t pos[3]={0.5,0.5,0.};
91	Double_t err[3]={3.,3.,5.};
92	initVertex = new AliESDVertex(pos,err);
93	}
94	invtx->Close();
95	delete invtx;
96	AliVertexerTracks *vertexer = new AliVertexerTracks();
97	vertexer->SetVtxStart(initVertex);
98	vertexer->SetDebug(0);
99	delete initVertex;
100	initVertex=0;
101
102	/*** The TREE for ESD is born ***/
103	TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
225a9242	104	AliESD event=0; AliESD eventTPCin=0;
970d08c3	105	esdTree->Branch("ESD","AliESD",&event);
	106
	107	if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
	108	if(lastEvent>rl->GetNumberOfEvents()) lastEvent=rl->GetNumberOfEvents()-1;
	109	cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
	110
	111	TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
	112	AliESDVertex *vertexSPD = new AliESDVertex();
	113	Char_t zver[100];
	114	Double_t vtx[3]={0,0,0};
	115	Double_t sigmavtx[3]={0.07,0.07,0.1};
	116
	117
	118	//<---------------------------------- The Loop over events begins
	119	TStopwatch timer;
	120	Int_t trc;
	121	for(Int_t i=firstEvent; i<=lastEvent; i++) {
	122
225a9242	123	cout<<" Processing event number : "<<i<<endl;
970d08c3	124	AliESD *event = new AliESD();
	125	event->SetRunNumber(gAlice->GetRunNumber());
	126	event->SetEventNumber(i);
	127	event->SetMagneticField(gAlice->Field()->SolenoidField());
	128	rl->GetEvent(i);
	129
	130	//***** Primary vertex from SPD from file
	131	sprintf(zver,"Event%d/Vertex",i);
	132	vertexSPD = (AliESDVertex*)ppZ->Get(zver);
	133	if(!vertexSPD) {
	134	esdTree->Fill(); delete event;
	135	continue;
	136	}
	137	event->SetVertex(vertexSPD);
	138	vertexSPD->GetXYZ(vtx);
	139	vertexSPD->GetSigmaXYZ(sigmavtx);
	140
	141	//***** TPC tracking
	142	if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
	143	printf("exiting TPC tracker with code %d in event %d\n",trc,i);
	144	esdTree->Fill(); delete event;
	145	continue;
	146	}
a654e257	147
225a9242	148	// make a copy of the ESD at this stage
	149	eventTPCin = event;
	150
970d08c3	151	//***** ITS tracking
	152	itsTracker.AliTracker::SetVertex(vtx,sigmavtx);
	153	TTree *itsTree=itsl->TreeR();
	154	if (!itsTree) {
	155	cerr<<"Can't get the ITS cluster tree !\n";
	156	esdTree->Fill(); delete event;
a654e257	157	return;
970d08c3	158	}
	159	itsTracker.UnloadClusters();
	160	itsTracker.LoadClusters(itsTree);
	161	if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
	162	printf("exiting ITS tracker with code %d in event %d\n",trc,i);
	163	esdTree->Fill(); delete event;
	164	continue;
	165	}
225a9242	166
	167	// Bring kTPCin-tracks back to the TPC inner wall
	168	BackToTPCInnerWall(event,eventTPCin);
	169
	170	// refit inward in ITS:
	171	// - refit without vertex constraint
	172	// - propagate through beam pipe to local x = 0
	173	itsTracker.RefitInward(event);
970d08c3	174
970d08c3	175
	176	//***** Vertex from ESD tracks
	177	if(collcode==1) { // pp
	178	AliESDVertex *vertexTrks =
	179	(AliESDVertex*)vertexer->FindPrimaryVertex(event);
	180	event->SetPrimaryVertex(vertexTrks);
	181	}
	182
	183	esdTree->Fill();
	184	delete event;
	185
	186	}//<-----------------------------------The Loop over events ends here
	187	timer.Stop(); timer.Print();
	188
	189	// The AliESDs.root is born
	190	TFile *ef = TFile::Open("AliESDs.root","RECREATE");
	191	if (!ef \|\| !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
	192
225a9242	193	//Write the tree and close everything
970d08c3	194	esdTree->Write();
	195	delete esdTree;
	196	ef->Close();
	197
	198	delete vertexer;
	199	delete rl;
	200
	201	return;
056891c0	202	}
4201c420	203	//--------------------------------------------------------------------------
225a9242	204	void BackToTPCInnerWall(AliESD event,AliESD eventTPC) {
	205
	206	Int_t ntracks = eventTPC->GetNumberOfTracks();
	207	AliESDtrack *esdTrackTPC = 0;
4201c420	208
225a9242	209	// create relation between event and eventTPC
	210	Int_t labelsTPC[100000000];
	211	for(Int_t tr = 0; tr<ntracks; tr++) {
	212	esdTrackTPC = (AliESDtrack*)event->GetTrack(tr);
	213	labelsTPC[TMath::Abs(esdTrackTPC->GetLabel())] = tr;
	214	}
	215
	216	ntracks = event->GetNumberOfTracks();
4201c420	217	AliESDtrack *esdTrack = 0;
225a9242	218	esdTrackTPC = 0;
225a9242	219	Int_t indexTPC;
4201c420	220
4201c420	221	// loop on tracks
225a9242	222	for(tr = 0; tr<ntracks; tr++) {
	223	esdTrack = (AliESDtrack*)event->GetTrack(tr);
	224	// set to kITSout the tracks that don't have kTPCin
	225	// (they've been found by AliITStrackerSA)
	226	if(!(esdTrack->GetStatus()&AliESDtrack::kTPCin)) {
	227	esdTrack->SetStatus(AliESDtrack::kITSout);
	228	continue;
	229	}
	230
	231	// skip tracks that don't have kITSin
4201c420	232	if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) continue;
225a9242	233
	234	indexTPC = labelsTPC[TMath::Abs(esdTrack->GetLabel())];
	235	esdTrackTPC = (AliESDtrack*)eventTPC->GetTrack(indexTPC);
	236
	237	AliITStrackMI *itsTrack = 0;
4201c420	238	try {
225a9242	239	itsTrack = new AliITStrackMI(*esdTrackTPC);
4201c420	240	esdTrack = 0;
	241	}
	242	catch (const Char_t *msg) {
225a9242	243	Warning("ToTPCInnerWall",msg);
4201c420	244	continue;
4201c420	245	}
225a9242	246	itsTrack->UpdateESDtrack(AliESDtrack::kITSout);
4201c420	247	esdTrack = new AliESDtrack(*(itsTrack->GetESDtrack()));
225a9242	248
4201c420	249	delete itsTrack;
	250	} // end loop on tracks
	251
970d08c3	252	return;
4201c420	253	}
225a9242	254	//--------------------------------------------------------------------------