first checkin of Claus' code with a few modifications
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / testAnalysis.C
1 //
2 // Script to test the dN/dEta analysis using the dNdEtaAnalysis and
3 // dNdEtaCorrection classes. Note that there is a cut on the events,
4 // so the measurement will be biassed.
5 //
6
7 testAnalysis(Char_t* dataDir, Int_t nRuns=20) {
8
9   Char_t str[256];
10
11   gSystem->Load("../esdTrackCuts/libESDtrackQuality.so");
12   gSystem->Load("libdNdEta.so");
13
14   // ########################################################
15   // selection of esd tracks
16   ESDtrackQualityCuts* esdTrackCuts = new ESDtrackQualityCuts();    
17   esdTrackCuts->DefineHistograms(1);
18   
19   esdTrackCuts->SetMinNClustersTPC(50);
20   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
21   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
22   esdTrackCuts->SetRequireTPCRefit(kTRUE);
23   
24   esdTrackCuts->SetMinNsigmaToVertex(3);
25   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
26
27   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
28
29   // ########################################################
30   // definition of dNdEta objects
31
32   dNdEtaCorrection* dNdEtaMap = new dNdEtaCorrection();
33
34   dNdEtaMap->LoadHistograms("correction_map.root","dndeta_correction");
35   dNdEtaMap->RemoveEdges(2,0,2);
36
37   dNdEtaAnalysis* analyse = new dNdEtaAnalysis("dndeta");
38
39   // ########################################################
40   // get the data dir  
41   Char_t execDir[256];
42   sprintf(execDir,gSystem->pwd());
43   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
44   TList* dirList            = baseDir->GetListOfFiles();
45   Int_t nDirs               = dirList->GetEntries();
46   // go back to the dir where this script is executed
47   gSystem->cd(execDir);
48
49
50   // ########################################################
51   // definition of used pointers
52   TFile* esdFile;
53   TTree* esdTree;
54   TBranch* esdBranch;
55
56   AliESD* esd =0;
57
58   // ########################################################
59   // loop over runs
60   Int_t nRunCounter = 0;
61   for (Int_t r=1; r<=nDirs; r++) {
62
63     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
64     if (!presentDir->IsDirectory())
65       continue;
66     // check that the files are there
67     sprintf(str,"%s/%s",dataDir, presentDir->GetName());
68     if ((!gSystem->Which(str,"galice.root")) ||
69         (!gSystem->Which(str,"AliESDs.root"))) 
70       continue;
71     
72     if (nRunCounter++ >= nRuns)
73       break;    
74     
75     cout << "run #" << nRunCounter << endl;
76     
77     // #########################################################
78     // setup galice and runloader
79     if (gAlice) {
80       delete gAlice->GetRunLoader();
81       delete gAlice;
82       gAlice=0;
83     }
84
85     sprintf(str,"%s/run%d/galice.root",dataDir,r);
86     AliRunLoader* runLoader = AliRunLoader::Open(str);
87
88     runLoader->LoadgAlice();
89     gAlice = runLoader->GetAliRun();
90     runLoader->LoadKinematics();
91     runLoader->LoadHeader();
92
93     // #########################################################
94     // open esd file and get the tree
95
96     sprintf(str,"%s/run%d/AliESDs.root",dataDir,r);
97     // close it first to avoid memory leak
98     if (esdFile)
99       if (esdFile->IsOpen())
100         esdFile->Close();
101
102     esdFile = TFile::Open(str);
103     esdTree = (TTree*)esdFile->Get("esdTree");
104     if (!esdTree)
105       continue;
106     esdBranch = esdTree->GetBranch("ESD");
107     esdBranch->SetAddress(&esd);
108     if (!esdBranch)
109       continue;
110
111     // ########################################################
112     // Magnetic field
113     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
114     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
115
116     // ########################################################
117     // getting number of events
118     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
119     Int_t nESDEvents = esdBranch->GetEntries();
120     
121     if (nEvents!=nESDEvents) {
122       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
123       return;
124     }
125     // ########################################################
126     // loop over number of events
127     cout << " looping over events..." << endl;
128     for(Int_t i=1; i<nEvents; i++) {
129
130       esdBranch->GetEntry(i);
131       runLoader->GetEvent(i);            
132   
133
134       // ########################################################
135       // get the EDS vertex
136       AliESDVertex* vtxESD = esd->GetVertex();
137
138       Double_t vtx[3];
139       Double_t vtx_res[3];
140       vtxESD->GetXYZ(vtx);          
141     
142       vtx_res[0] = vtxESD->GetXRes();
143       vtx_res[1] = vtxESD->GetYRes();
144       vtx_res[2] = vtxESD->GetZRes();
145
146       // we need a good vertex 
147       //  => there will be a bias on the measurement, since this cuts away some events
148       if (strcmp(vtxESD->GetName(),"default")==0) 
149         continue;
150       if (vtx_res[2]==0 || vtx_res[2]>0.1)
151         continue;
152
153       // ########################################################
154       // loop over esd tracks      
155       Int_t nTracks = esd->GetNumberOfTracks();            
156       for (Int_t t=0; t<nTracks; t++) {
157         AliESDtrack* esdTrack = esd->GetTrack(t);      
158         
159         if (!esdTrackCuts->AcceptTrack(esdTrack))
160           continue;
161         
162         AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
163         
164         if (tpcTrack->GetAlpha()==0) {
165           cout << " Warning esd track alpha = 0" << endl;
166           continue;
167         }
168
169         Float_t theta = tpcTrack->Theta();
170         Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
171         Float_t correction = dNdEtaMap->GetCorrection(vtx[2], eta);
172         
173         dNdEtaMap->FillMeas(vtx[2], eta);
174
175         analyse   ->FillTrack(vtx[2], eta, correction);
176         
177       } // end of track loop
178       analyse->FillEvent(vtx[2]);
179
180     } // end  of event loop
181   } // end of run loop
182
183   analyse->Finish();
184
185   TFile* fout = new TFile("out.root","RECREATE");
186   
187   esdTrackCuts->SaveHistograms("esd_tracks_cuts");
188   dNdEtaMap->SaveHistograms();
189   analyse   ->SaveHistograms();
190   
191   fout->Write();
192   fout->Close();
193
194 }