Corrected cast
[u/mrichter/AliRoot.git] / PWG0 / ptSpectra / makeCorrectionPtEtaVSz.C
1 /* $Id$ */
2
3 //
4 // Script to calculate PtvsEta correction map at different z of vtx.
5 // 
6
7 makeCorrectionPtEtaVSz(Char_t* dataDir, Int_t nRuns=20) {
8
9   Char_t str[256];
10 static const Int_t  NZbin=5;
11 Float_t Zbin[]={-10., -3.,0., 3., 10.};
12
13   gSystem->Load("../libPWG0base.so");
14   // ########################################################
15   // selection of esd tracks
16   AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts();    
17   esdTrackCuts->DefineHistograms(1);
18   esdTrackCuts->SetMinNClustersTPC(50);
19   esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
20   esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
21   esdTrackCuts->SetRequireTPCRefit(kTRUE);
22   esdTrackCuts->SetMinNsigmaToVertex(3);
23   esdTrackCuts->SetAcceptKingDaughters(kFALSE);
24   AliLog::SetClassDebugLevel("AliESDtrackCuts",1);
25   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
26
27   // ########################################################
28   // definition of PtEta correction objects
29
30
31 CorrectionMatrix2D *PtEtaMap[NZbin];
32 for(Int_t i=0; i<NZbin-1; i++)
33   PtEtaMap[i] = new CorrectionMatrix2D(Form("PtvsEta%d",i),Form("%f<z_{vtx}<%f",Zbin[i],Zbin[i+1]),80,0.,10.,120,-2.,2.);
34
35   PtEtaMap[NZbin-1] = new CorrectionMatrix2D("xxx","xxx",80,0.,10.,120,-2.,2.);
36
37
38   // ########################################################
39   // get the data dir  
40   Char_t execDir[256];
41   sprintf(execDir,gSystem->pwd());
42   TSystemDirectory* baseDir = new TSystemDirectory(".",dataDir);
43   TList* dirList            = baseDir->GetListOfFiles();
44   Int_t nDirs               = dirList->GetEntries();
45   // go back to the dir where this script is executed
46   gSystem->cd(execDir);
47
48   // ########################################################
49   // definition of used pointers
50   TFile* esdFile;
51   TTree* esdTree;
52   TBranch* esdBranch;
53
54   AliESD* esd =0;
55   // ########################################################
56   // loop over runs
57   Int_t nRunCounter = 0;
58   for (Int_t r=0; r<nDirs; r++) {
59
60     TSystemFile* presentDir = (TSystemFile*)dirList->At(r);
61     if (!presentDir || !presentDir->IsDirectory())
62       continue;
63     // check that the files are there
64     TString currentDataDir;
65     currentDataDir.Form("%s/%s",dataDir, presentDir->GetName());
66     cout << "Processing directory " << currentDataDir.Data() << endl;
67     if ((!gSystem->Which(currentDataDir,"galice.root")) ||
68           (!gSystem->Which(currentDataDir,"AliESDs.root")))
69       continue;
70
71     if (nRunCounter++ >= nRuns)
72       break;
73
74     cout << "run #" << nRunCounter << endl;
75
76     // #########################################################
77     // setup galice and runloader
78     if (gAlice) {
79       delete gAlice->GetRunLoader();
80       delete gAlice;
81       gAlice=0;
82     }
83
84     sprintf(str,"%s/galice.root",currentDataDir.Data());
85     AliRunLoader* runLoader = AliRunLoader::Open(str);
86
87     runLoader->LoadgAlice();
88     gAlice = runLoader->GetAliRun();
89     runLoader->LoadKinematics();
90     runLoader->LoadHeader();
91
92    
93     // #########################################################
94     // open esd file and get the tree
95
96     sprintf(str,"%s/AliESDs.root",currentDataDir.Data());
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=0; i<nEvents; i++) {
129
130       esdBranch->GetEntry(i);
131       runLoader->GetEvent(i);              
132
133       // ########################################################
134       // get the EDS vertex
135       AliESDVertex* vtxESD = esd->GetVertex();
136
137       Double_t vtx[3];
138       Double_t vtx_res[3];
139       vtxESD->GetXYZ(vtx);          
140
141     
142       vtx_res[0] = vtxESD->GetXRes();
143       vtx_res[1] = vtxESD->GetYRes();
144       vtx_res[2] = vtxESD->GetZRes();
145
146       // the vertex should be reconstructed
147       if (strcmp(vtxESD->GetName(),"default")==0) 
148         continue;
149
150       // the resolution should be reasonable???
151       if (vtx_res[2]==0 || vtx_res[2]>0.1)
152         continue;
153
154 Int_t IZbin=NZbin-1;
155           for(Int_t l=0; l<NZbin-1; l++)
156             if((vtx[3]<Zbin[l+1])&&((vtx[3]>Zbin[l])))IZbin=l;
157
158
159
160       // ########################################################
161       // loop over mc particles
162       AliStack* particleStack = runLoader->Stack();
163       Int_t nPrim    = particleStack->GetNprimary();
164
165       for (Int_t i_mc=0; i_mc<nPrim; i_mc++) {
166       
167         TParticle* part = particleStack->Particle(i_mc);      
168         if (!part || strcmp(part->GetName(),"XXX")==0) 
169           continue;
170       
171         TParticlePDG* pdgPart = part->GetPDG();
172
173         Bool_t prim = kFALSE;
174         // check if it's a primary particle - is there a better way ???
175         if ((part->GetFirstDaughter() >= nPrim) || (part->GetFirstDaughter()==-1)) {
176           if (TMath::Abs(pdgPart->PdgCode())>10 && pdgPart->PdgCode()!=21 && strcmp(pdgPart->ParticleClass(),"Unknown")!=0)
177             prim = kTRUE;
178         }
179         if (!prim)
180           continue;
181
182         if (pdgPart->Charge()==0)
183           continue;
184
185         if (prim)
186             PtEtaMap[IZbin]->FillGene(part->Pt(), part->Eta()); 
187         
188       }// end of mc particle
189
190       // ########################################################
191       // loop over esd tracks      
192       Int_t nTracks = esd->GetNumberOfTracks();            
193       for (Int_t t=0; t<nTracks; t++) {
194         AliESDtrack* esdTrack = esd->GetTrack(t);      
195         
196         // cut the esd track?
197         if (!esdTrackCuts->AcceptTrack(esdTrack))
198           continue;
199
200         AliTPCtrack* tpcTrack = new AliTPCtrack(*esdTrack);
201         if (tpcTrack->GetAlpha()==0) {
202           cout << " Warning esd track alpha = 0" << endl;
203           continue;
204         }
205
206         Float_t theta = tpcTrack->Theta();
207         Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
208         Float_t pt=  = tpcTrack->Pt();
209
210
211         PtEtaMap[IZbin]->FillMeas(eta, pt);     
212
213       } // end of track loop
214     } // end  of event loop
215   } // end of run loop
216
217   
218   for(Int_t l=0; l<NZbin-1; l++)
219   {
220   PtEtaMap[l]->Divide();  
221   PtEtaMap[l]->RemoveEdges();  
222   }
223   TFile* fout = new TFile("PtEtaCorrectionMap.root","RECREATE");
224   
225   esdTrackCuts->SaveHistograms("EsdTrackCuts");
226   for(Int_t l=0; l<NZbin-1; l++) PtEtaMap[l]->SaveHistograms();
227   
228   fout->Write();
229   fout->Close();
230
231 }