o) Added DrawHistograms function
[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     TString currentDataDir;
68     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
69     if ((!gSystem->Which(currentDataDir.Data(),"galice.root")) ||
70         (!gSystem->Which(currentDataDir.Data(),"AliESDs.root")))
71       continue;
72     
73     if (nRunCounter++ >= nRuns)
74       break;    
75     
76     cout << "run #" << nRunCounter << endl;
77     
78     // #########################################################
79     // setup galice and runloader
80     if (gAlice) {
81       delete gAlice->GetRunLoader();
82       delete gAlice;
83       gAlice=0;
84     }
85
86     sprintf(str,"%s/galice.root",currentDataDir.Data());
87     AliRunLoader* runLoader = AliRunLoader::Open(str);
88
89     runLoader->LoadgAlice();
90     gAlice = runLoader->GetAliRun();
91     runLoader->LoadKinematics();
92     runLoader->LoadHeader();
93
94     // #########################################################
95     // open esd file and get the tree
96
97     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
98     // close it first to avoid memory leak
99     if (esdFile)
100       if (esdFile->IsOpen())
101         esdFile->Close();
102
103     esdFile = TFile::Open(str);
104     esdTree = (TTree*)esdFile->Get("esdTree");
105     if (!esdTree)
106       continue;
107     esdBranch = esdTree->GetBranch("ESD");
108     esdBranch->SetAddress(&esd);
109     if (!esdBranch)
110       continue;
111
112     // ########################################################
113     // Magnetic field
114     AliTracker::SetFieldMap(gAlice->Field(),kTRUE); // kTRUE means uniform magnetic field
115     //AliKalmanTrack::SetConvConst(1000/0.299792458/5.);
116
117     // ########################################################
118     // getting number of events
119     Int_t nEvents = (Int_t)runLoader->GetNumberOfEvents();
120     Int_t nESDEvents = esdBranch->GetEntries();
121     
122     if (nEvents!=nESDEvents) {
123       cout << " Different number of events from runloader and esdtree!!!" << nEvents << " / " << nESDEvents << endl;
124       return;
125     }
126     // ########################################################
127     // loop over number of events
128     cout << " looping over events..." << endl;
129     for(Int_t i=1; i<nEvents; i++) {
130
131       esdBranch->GetEntry(i);
132       runLoader->GetEvent(i);            
133   
134
135       // ########################################################
136       // get the EDS vertex
137       AliESDVertex* vtxESD = esd->GetVertex();
138
139       Double_t vtx[3];
140       Double_t vtx_res[3];
141       vtxESD->GetXYZ(vtx);          
142     
143       vtx_res[0] = vtxESD->GetXRes();
144       vtx_res[1] = vtxESD->GetYRes();
145       vtx_res[2] = vtxESD->GetZRes();
146
147       // we need a good vertex 
148       //  => there will be a bias on the measurement, since this cuts away some events
149       if (strcmp(vtxESD->GetName(),"default")==0) 
150         continue;
151       if (vtx_res[2]==0 || vtx_res[2]>0.1)
152         continue;
153
154       // ########################################################
155       // loop over esd tracks      
156       Int_t nTracks = esd->GetNumberOfTracks();            
157       for (Int_t t=0; t<nTracks; t++) {
158         AliESDtrack* esdTrack = esd->GetTrack(t);      
159         
160         if (!esdTrackCuts->AcceptTrack(esdTrack))
161           continue;
162         
163         AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
164         
165         if (tpcTrack->GetAlpha()==0) {
166           cout << " Warning esd track alpha = 0" << endl;
167           continue;
168         }
169
170         Float_t theta = tpcTrack->Theta();
171         Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
172         Float_t correction = dNdEtaMap->GetCorrection(vtx[2], eta);
173         
174         dNdEtaMap->FillMeas(vtx[2], eta);
175
176         analyse   ->FillTrack(vtx[2], eta, correction);
177         
178       } // end of track loop
179       analyse->FillEvent(vtx[2]);
180
181     } // end  of event loop
182   } // end of run loop
183
184   analyse->Finish();
185
186   TFile* fout = new TFile("out.root","RECREATE");
187   
188   esdTrackCuts->SaveHistograms("esd_tracks_cuts");
189   dNdEtaMap->SaveHistograms();
190   analyse   ->SaveHistograms();
191   
192   fout->Write();
193   fout->Close();
194
195 }