Updated version of the fast TPC simulation which works with the ESD IO (A.Dainese)
[u/mrichter/AliRoot.git] / TPC / AliBarrelRec_TPCparam.C
1 void AliBarrelRec_TPCparam(Int_t firstEvent=0,Int_t lastEvent=0) {
2   //
3   // Macro to create AliESDs.root using parametrized TPC tracking
4   // and AliITStrackerV2
5   //
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)
12   //
13   // A. Dainese - LNL
14   //
15
16
17    /**** Initialization of the NewIO *******/
18
19    if (gAlice) {
20       delete gAlice->GetRunLoader();
21       delete gAlice; 
22       gAlice=0;
23    }
24
25    AliRunLoader *rl = AliRunLoader::Open("galice.root");
26    if (rl == 0x0) {
27       cerr<<"Can not open session"<<endl;
28       return;
29    }
30    Int_t retval = rl->LoadgAlice();
31    if (retval) {
32       cerr<<"LoadgAlice returned error"<<endl;
33       delete rl;
34       return;
35    }
36    retval = rl->LoadHeader();
37    if (retval) {
38       cerr<<"LoadHeader returned error"<<endl;
39       delete rl;
40       return;
41    }
42    gAlice=rl->GetAliRun();
43        
44
45    TDatabasePDG *DataBase = TDatabasePDG::Instance();
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    Int_t collcode = 1; // pp collisions
56    AliTPCtrackerParam tpcTrackerPar(collcode,fieval);
57    tpcTrackerPar.Init();
58
59    //**** Switch on the PID class (mandatory!)
60    AliPID pid;
61
62    /**** The ITS corner ********************/
63
64    AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
65    if (itsl == 0x0) {
66       cerr<<"Cannot get the ITS loader"<<endl;
67       return;
68    }
69    itsl->LoadRecPoints("read");
70
71    AliITS *dITS = (AliITS*)gAlice->GetDetector("ITS");
72    if (!dITS) {
73       cerr<<"Cannot find the ITS detector !"<<endl;
74       return;
75    }
76    AliITSgeom *geom = dITS->GetITSgeom();
77
78    //An instance of the ITS tracker
79    AliITStrackerV2 itsTracker(geom);
80    Int_t ITSclusters[6] = {1,1,1,1,1,1};
81    itsTracker.SetLayersNotToSkip(ITSclusters);
82
83    /***** The TREE is born *****/
84    
85    TTree *esdTree=new TTree("esdTree","Tree with ESD objects");
86    AliESD *event=0;
87    esdTree->Branch("ESD","AliESD",&event);
88    
89    if(firstEvent>rl->GetNumberOfEvents()) firstEvent=rl->GetNumberOfEvents()-1;
90    if(lastEvent>rl->GetNumberOfEvents())  lastEvent=rl->GetNumberOfEvents()-1;
91    cout<<" Number of events: "<<1+lastEvent-firstEvent<<endl;
92    
93    TFile *ppZ = TFile::Open("ITS.Vertex.root"); // z vertices from SPD
94    AliESDVertex *myvertex = new AliESDVertex();
95    Char_t zver[100];
96    Double_t vtx[3];
97    Double_t cvtx[6];
98
99
100    //<----------------------------------The Loop over events begins
101    TStopwatch timer;
102    Int_t trc;
103    for(Int_t i=firstEvent; i<=lastEvent; i++) { 
104      
105      cerr<<" Processing event number : "<<i<<endl;
106      AliESD *event = new AliESD(); 
107      event->SetRunNumber(gAlice->GetRunNumber());
108      event->SetEventNumber(i);
109      event->SetMagneticField(gAlice->Field()->SolenoidField());
110      rl->GetEvent(i);
111
112     //***** Primary vertex 
113      sprintf(zver,"Event%d/Vertex",i);
114      myvertex = (AliESDVertex*)ppZ->Get(zver);
115      if(!myvertex) {
116        esdTree->Fill(); delete event;
117        continue;
118      }      
119      event->SetVertex(myvertex);
120      myvertex->GetXYZ(vtx);
121      myvertex->GetCovMatrix(cvtx);
122
123      //***** TPC tracking
124      if ( (trc=tpcTrackerPar.BuildTPCtracks(event)) ) {
125        printf("exiting TPC tracker with code %d in event %d\n",trc,i);
126        esdTree->Fill(); delete event;
127        continue;
128      }
129
130      //***** ITS tracking
131      itsTracker.SetVertex(vtx,cvtx);
132      TTree *itsTree=itsl->TreeR();
133      if (!itsTree) {
134         cerr<<"Can't get the ITS cluster tree !\n";
135         esdTree->Fill(); delete event;
136         return;
137      }     
138      itsTracker.UnloadClusters();
139      itsTracker.LoadClusters(itsTree);
140      if ( (trc=itsTracker.Clusters2Tracks(event)) ) {
141        printf("exiting ITS tracker with code %d in event %d\n",trc,i);
142        esdTree->Fill(); delete event;
143        continue;
144      }
145
146
147      esdTree->Fill();
148      delete event;
149
150    }//<-----------------------------------The Loop over events ends here
151    timer.Stop(); timer.Print();
152
153    //        The AliESDs.root is born
154    TFile *ef = TFile::Open("AliESDs.root","RECREATE"); 
155    if (!ef || !ef->IsOpen()) {cerr<<"Can't open AliESDs.root !\n"; return;}
156
157    esdTree->Write(); //Write the TREE and close everything
158    delete esdTree;
159    ef->Close();
160
161    delete rl;
162
163    return;
164 }