changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TOF / pp7 / TOFSpectrappAnalysis.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ///////////////////////////////////////////////////////////////////////////
17 //                                                                       //
18 //                                                                       //
19 // Analysis for identified charged hadron spectra. TOF                   //
20 //                                                                       //
21 //                                                                       //
22 ///////////////////////////////////////////////////////////////////////////
23
24 #include "Riostream.h"
25 #include "TChain.h"
26 #include "TTree.h"
27 #include "TH1F.h"
28 #include "TH2F.h"
29 #include "TH3F.h"
30 #include "TList.h"
31 #include "TMath.h"
32 #include "TCanvas.h"
33 #include "TObjArray.h"
34 #include "TF1.h"
35 #include "TFile.h"
36
37 #include "AliAnalysisTaskSE.h"
38 #include "AliAnalysisManager.h"
39
40 #include "AliHeader.h"
41 #include "AliGenPythiaEventHeader.h"
42 #include "AliGenCocktailEventHeader.h"
43
44 #include "AliPID.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliESDVertex.h"
47 #include "AliESDEvent.h"
48 #include "AliESDInputHandler.h"
49 #include "AliESDtrack.h"
50 #include "AliESDpid.h"
51
52 #include "AliMCEventHandler.h"
53 #include "AliMCEvent.h"
54 #include "AliStack.h"
55
56 #include "AliAnalysisFilter.h"
57 #include "AliCFContainer.h"
58 #include "AliCFFrame.h"
59 #include "AliCFGridSparse.h"
60 #include "TDatabasePDG.h"
61 #include "TROOT.h"
62 #include "TSystem.h"
63 #include "AliMultiplicity.h"
64
65 #include "AliLog.h"
66
67 #include "TOFSpectrappAnalysis.h"
68 //cambio
69 #include "AliESDUtils.h"
70 //fine cambio
71 ClassImp(TOFSpectrappAnalysis)
72
73 //________________________________________________________________________
74 TOFSpectrappAnalysis::TOFSpectrappAnalysis() 
75 : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fESDtrackCuts(0), fESDpid(0x0),
76   fMCtrue(0),
77   fAlephParameters(),
78   spdCorr(-1),
79   multiplicity(0),
80   ZPrimVertex(0),
81   frun(0),
82   frunOld(0),
83   fLoadOCDB(kTRUE),
84   correctTExp(kTRUE),
85   calibrateESD(kTRUE), useT0TOF(kTRUE), timeResolution(96.), tuneTOFMC(kFALSE),
86   fTreeTrack(0),
87   //fTreeEv(0),
88   hNumEv(0),
89   tofCalib(0),
90   t0maker(0),
91   fDCAXY(-999), fDCAZ(-999), kselimppar(-999), fmatch(-999), fmom(-999), flength(-999), fsign(-999), ftimetof(-999), fexptimeel(-999), fexptimemu(-999), fexptimepi(-999), fexptimeka(-999), fexptimepr(-999), ftofchan(-999), feta(-999), fphi(-999), fmomtrasv(-999), t0track(-999), t0trackSigma(-999), sigmael(-999), sigmamu(-999), sigmapi(-999), sigmaka(-999), sigmapr(-999), TPCSignal(-999), TPCSigmaPI(-999), TPCSigmaKA(-999), TPCSigmaPR(-999), TOFlabel0(-999), TOFlabel1(-999), TOFlabel2(-999)    
92 {
93   // default Constructor
94   for(Int_t i=0; i<5;i++){r1[i]=-999;}
95 }
96
97
98 //________________________________________________________________________
99 TOFSpectrappAnalysis::TOFSpectrappAnalysis(const char *name) 
100   : AliAnalysisTaskSE(name), fESD(0),fESDtrackCuts(0),fESDpid(0x0),
101     fMCtrue(0),
102     fAlephParameters(),
103     spdCorr(-1),
104     multiplicity(0),
105     ZPrimVertex(-999),
106     frun(0),
107     frunOld(0),
108     fLoadOCDB(kTRUE),
109     correctTExp(kTRUE),
110     calibrateESD(kTRUE), useT0TOF(kTRUE), timeResolution(96.), tuneTOFMC(kFALSE),
111     fTreeTrack(0),
112     //fTreeEv(0),
113     hNumEv(0),
114     tofCalib(0),
115     t0maker(0),
116     fDCAXY(-999), fDCAZ(-999), kselimppar(-999), fmatch(-999), fmom(-999), flength(-999), fsign(-999), ftimetof(-999), fexptimeel(-999), fexptimemu(-999), fexptimepi(-999), fexptimeka(-999), fexptimepr(-999), ftofchan(-999), feta(-999), fphi(-999), fmomtrasv(-999), t0track(-999), t0trackSigma(-999), sigmael(-999), sigmamu(-999), sigmapi(-999), sigmaka(-999), sigmapr(-999), TPCSignal(-999), TPCSigmaPI(-999), TPCSigmaKA(-999), TPCSigmaPR(-999), TOFlabel0(-999), TOFlabel1(-999), TOFlabel2(-999)    
117 {
118   
119   //
120   // standard constructur which should be used
121   //
122   Printf("*** CONSTRUCTOR CALLED ****");
123   // 
124   /* real */
125   for(Int_t i=0; i<5;i++){r1[i]=-999;}
126
127
128   fAlephParameters[0] = 0.0283086;
129   fAlephParameters[1] = 2.63394e+01;
130   fAlephParameters[2] = 5.04114e-11;
131   fAlephParameters[3] = 2.12543e+00;
132   fAlephParameters[4] = 4.88663e+00;
133   //
134   // initialize PID object
135   //
136
137   
138   //tofCalib = new AliTOFcalib();
139
140   //cambio
141
142   //fESDpid = new AliESDpid();
143
144   //fine cambio
145   //t0maker = new AliTOFT0maker(fESDpid, tofCalib); 
146   //
147   // create track cuts
148   //
149   AliESDtrackCuts* ESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
150   fESDtrackCuts = (AliESDtrackCuts*) ESDtrackCuts->GetStandardITSTPCTrackCuts2010(kTRUE);
151
152   //fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
153
154   //
155   //Initialize();
156   // Output slot #0 writes into a TList container
157   //DefineOutput(1, TTree::Class());
158   DefineOutput(1, TTree::Class());
159   DefineOutput(2, TH1D::Class());
160
161 }
162
163
164 //________________________________________________________________________
165 Int_t TOFSpectrappAnalysis::Mult()
166 {
167   AliESDtrackCuts* esdTrackCutsMult = new AliESDtrackCuts;
168   
169   // TPC
170   esdTrackCutsMult->SetMinNClustersTPC(70);//ok
171   esdTrackCutsMult->SetMaxChi2PerClusterTPC(4);//ok
172   esdTrackCutsMult->SetAcceptKinkDaughters(kFALSE);//ok
173   esdTrackCutsMult->SetRequireTPCRefit(kTRUE);//ok
174   // ITS
175   esdTrackCutsMult->SetRequireITSRefit(kTRUE);//ok
176   esdTrackCutsMult->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
177                                              AliESDtrackCuts::kAny);//ok
178   
179   esdTrackCutsMult->SetMaxDCAToVertexXYPtDep("0.0182+0.0350/pt^1.01");
180   
181   esdTrackCutsMult->SetMaxDCAToVertexZ(2);//ok
182   
183   esdTrackCutsMult->SetDCAToVertex2D(kFALSE);//ok
184   esdTrackCutsMult->SetRequireSigmaToVertex(kFALSE);//ok
185   
186   esdTrackCutsMult->SetEtaRange(-0.8,+0.8);
187   esdTrackCutsMult->SetPtRange(0.15, 1e10);
188   
189   // Here is what we finally want:
190   Int_t multipli = esdTrackCutsMult->CountAcceptedTracks(fESD);
191   delete esdTrackCutsMult;
192   return multipli;
193   
194  }
195
196
197 //________________________________________________________________________
198 void TOFSpectrappAnalysis::UserCreateOutputObjects() 
199 {
200   // Create histograms
201   // Called once
202   
203   OpenFile(1);
204
205   fTreeTrack = new TTree("TreeTrack","Track Properties");
206   fTreeTrack->Branch("spdCorr",&spdCorr,"spdCorr/F");  
207   fTreeTrack->Branch("multiplicity",&multiplicity,"multiplicity/I");
208   fTreeTrack->Branch("DCAXY",&fDCAXY,"fDCAXY/F");  
209   fTreeTrack->Branch("DCAZ",&fDCAZ,"fDCAZ/F");  
210   //fTreeTrack->Branch("ZPrimVertex",&ZPrimVertex,"ZPrimVertex/D");  
211   //fTreeTrack->Branch("kselimppar",&kselimppar,"kselimppar/I");  
212   fTreeTrack->Branch("fmatch",&fmatch,"fmatch/I");
213   fTreeTrack->Branch("fmom",&fmom,"fmom/D");
214   fTreeTrack->Branch("length",&flength,"length/D");
215   fTreeTrack->Branch("sign",&fsign,"sign/I");
216   fTreeTrack->Branch("timetof",&ftimetof,"timetof/D");
217   fTreeTrack->Branch("exptimeel",&fexptimeel,"exptimeel/D");
218   fTreeTrack->Branch("exptimemu",&fexptimemu,"exptimemu/D");
219   fTreeTrack->Branch("exptimepi",&fexptimepi,"exptimepi/D");
220   fTreeTrack->Branch("exptimeka",&fexptimeka,"exptimeka/D");
221   fTreeTrack->Branch("exptimepr",&fexptimepr,"exptimepr/D");
222   fTreeTrack->Branch("tofchan",&ftofchan,"tofchan/I"); 
223   fTreeTrack->Branch("eta",&feta,"eta/D");
224   fTreeTrack->Branch("phi",&fphi,"phi/D");
225   fTreeTrack->Branch("momtrasv",&fmomtrasv,"momtrasv/D");
226   fTreeTrack->Branch("t0trackSigma",&t0trackSigma,"t0trackSigma/F");
227   fTreeTrack->Branch("t0track",&t0track,"t0track/F");
228   fTreeTrack->Branch("sigmael",&sigmael,"sigmael/D");
229   fTreeTrack->Branch("sigmamu",&sigmamu,"sigmamu/D");
230   fTreeTrack->Branch("sigmapi",&sigmapi,"sigmapi/D");
231   fTreeTrack->Branch("sigmaka",&sigmaka,"sigmaka/D");
232   fTreeTrack->Branch("sigmapr",&sigmapr,"sigmapr/D");
233  //  //fTreeTrack->Branch("TPCsignal",&TPCSignal,"TPCsignal/D");
234 //   //fTreeTrack->Branch("TPCSigmaPI",&TPCSigmaPI,"TPCSigmaPI/F");
235 //   //fTreeTrack->Branch("TPCSigmaKA",&TPCSigmaKA,"TPCSigmaKA/F");
236 //   //fTreeTrack->Branch("TPCSigmaPR",&TPCSigmaPR,"TPCSigmaPR/F");
237 //   fTreeTrack->Branch("r10",&r1[0],"r10/D");
238 //   fTreeTrack->Branch("r11",&r1[1],"r11/D");
239 //   fTreeTrack->Branch("r12",&r1[2],"r12/D");
240 //   fTreeTrack->Branch("r13",&r1[3],"r13/D");
241 //   fTreeTrack->Branch("r14",&r1[4],"r14/D");
242 //   // fTreeTrack->Branch("TOFlabel0",&TOFlabel0,"TOFlabel0/I");
243 // //   fTreeTrack->Branch("TOFlabel1",&TOFlabel1,"TOFlabel1/I");
244 // //   fTreeTrack->Branch("TOFlabel2",&TOFlabel2,"TOFlabel2/I");
245
246  
247   PostData(1, fTreeTrack );
248
249   // fTreeEv = new TTree("TreeEv","Event Properties");
250 //   fTreeEv->Branch("multiplicity",&multiplicity,"multiplicity/I");  
251 //   fTreeEv->Branch("spdCorr",&spdCorr,"spdCorr/F");  
252 //   //fTreeEv->Branch("ZPrimVertex",&ZPrimVertex,"ZPrimVertex/D"); 
253   
254
255   OpenFile(2);
256   hNumEv=new TH1D("NemEv","NemEv",4,1,5);
257   PostData(2, hNumEv);
258
259 }
260
261 //________________________________________________________________________
262 void TOFSpectrappAnalysis::UserExec(Option_t *) 
263 {
264   //
265   // main event loop
266   //
267  
268   
269   multiplicity=-999, ZPrimVertex=-999, frun=-999, spdCorr=-1.0;
270
271   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
272   if (!fESD) {
273     Printf("ERROR: fESD not available");
274     return;
275   }
276   
277   if (!fESDtrackCuts) {
278     Printf("ERROR: fESDtrackCuts not available");
279     return;
280   }
281
282   hNumEv->Fill(1);
283
284   //cambio
285   AliESDInputHandler* esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
286   if (esdH)
287     fESDpid = esdH->GetESDpid();
288   //fine cambio
289   //
290   // check if event is selected by physics selection class
291   //
292   Bool_t isSelected = kFALSE;
293   isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected()& AliVEvent::kMB);
294
295   if (!isSelected) {
296     return;
297   }
298   
299   hNumEv->Fill(2);
300
301   //
302   // monitor vertex position
303   //
304   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks(); //! Primary vertex estimated using ESD tracks
305   if(vertex->GetNContributors()<1) { // # of tracklets/tracks used for the estimate
306     // SPD vertex
307     vertex = fESD->GetPrimaryVertexSPD(); //! Primary vertex estimated by the SPD
308     if(vertex->GetNContributors()<1) vertex = 0x0;
309   }  
310
311   if (!vertex) {
312     return;
313   } 
314   
315   hNumEv->Fill(3);
316
317   if (vertex) {ZPrimVertex=vertex->GetZ();}
318   if(TMath::Abs(ZPrimVertex)>10){return;}
319   
320   hNumEv->Fill(4);
321   
322   multiplicity=Mult();
323   
324   //multiplicity as defined by Marek
325   const AliMultiplicity *mult = fESD->GetMultiplicity();
326   Float_t nClusters[6]={0.0,0.0,0.0,0.0,0.0,0.0};
327   for(Int_t ilay=0; ilay<6; ilay++)
328     {
329       nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
330     }
331   //cambio
332   spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],vertex->GetZ());
333   //fine cambio
334   //fTreeEv->Fill();
335   
336
337   //TOF settings done in the TOF tender
338   
339   frun = fESD->GetRunNumber();
340  //  if(frun==frunOld){fLoadOCDB=kFALSE;}else {fLoadOCDB=kTRUE;}
341 //   //if (tuneTOFMC) calibrateESD = kFALSE;
342 //   Double_t *T0TOF;
343 //   if(fLoadOCDB){
344 //   AliCDBManager *cdb = AliCDBManager::Instance();
345 //   cdb->SetDefaultStorage("alien://folder=/alice/data/2010/OCDB");
346 //   //cdb->SetDefaultStorage("raw://");
347 //   cdb->SetRun(frun);}
348   
349 //   /* init TOF calibration */
350 //   if (correctTExp)
351 //     tofCalib->SetCorrectTExp(kTRUE);
352 //   tofCalib->Init();
353   
354 //   /* init TOF T0-maker */
355 //   t0maker->SetTimeResolution(timeResolution);
356   
357 //   /* calibrate ESD */
358 //   if (calibrateESD)
359 //     tofCalib->CalibrateESD(fESD);
360    
361 //   T0TOF=t0maker->ComputeT0TOF(fESD);// calcola il t0 solo col tof in 10 bin di pt e lo setta in AliTOFPidResponse 
362
363 //   //scrive i valori precedenti nel TOFHeader
364 //   t0maker->WriteInESD(fESD); 
365
366 //   //setta T0_TOF in AliTOFPidResponse ovvero i valori settati in AliTOFHeader. Se non ci sono setta T0spread
367 //   fESDpid->SetTOFResponse(fESD,AliESDpid::kTOF_T0); 
368
369
370 //   fESDpid->MakePID(fESD,kFALSE); //calcola la sigma e le gi in piĆ¹ definisce la flag kTOFmismatch
371   
372
373   //
374   // track loop
375   //
376   
377   //const Float_t kNsigmaCut = 3;
378   //const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
379   //Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
380   
381   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
382     
383     AliESDtrack *track =fESD->GetTrack(i); 
384     if (!track){continue;}
385     // start TOF analysis
386     
387     fDCAXY=-999, fDCAZ=-999, kselimppar=-999, fmatch=-999, fmom=-999, flength=-999, fsign=-999, ftimetof=-999, fexptimeel=-999, fexptimeel=-999, fexptimepi=-999, fexptimeka=-999, fexptimepr=-999, ftofchan=-999, feta=-999, fphi=-999, fmomtrasv=-999, sigmael=-999, sigmamu=-999, sigmapi=-999, t0track=-999, t0trackSigma=-999, sigmaka=-999, sigmapr=-999, TPCSignal=-999, TPCSigmaPI=-999, TPCSigmaKA=-999, TPCSigmaPR=-999, r1[0]=-999,r1[1]=-999,r1[2]=-999,r1[3]=-999,r1[4]=-999;
388     
389     if (!fESDtrackCuts->AcceptTrack(track)) {continue;}
390     if (SelectOnImpPar(track)) {kselimppar=1;}else {kselimppar=0;}
391     if ((track->GetStatus()&AliESDtrack::kTOFout)==0) {continue;}
392     if ((track->GetStatus()&AliESDtrack::kTIME)==0) {continue;}  
393     if (!(track->GetStatus()&AliESDtrack::kTOFmismatch)==0) {fmatch=0;}else {fmatch=1;}  
394     track->GetImpactParameters(fDCAXY, fDCAZ);
395     track->GetTOFpid(r1);
396     fmom=track->GetP();  // This function returns the track momentum
397     flength=track->GetIntegratedLength();
398     fsign=track->GetSign();
399     ftimetof=track->GetTOFsignal();
400     Double_t inttime[5]; 
401     track->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
402     fexptimeel=inttime[0];
403     fexptimemu=inttime[1];
404     fexptimepi=inttime[2];
405     fexptimeka=inttime[3];
406     fexptimepr=inttime[4];  
407     ftofchan=track->GetTOFCalChannel(); // Channel Index of the TOF Signal
408     feta=track->Eta();// return pseudorapidity return -TMath::Log(TMath::Tan(0.5 * Theta())); 
409     fphi=track->Phi()* 180 / TMath::Pi();// Returns the azimuthal angle of momentum  0 <= phi < 2*pi
410     fmomtrasv=track->Pt(); 
411     t0track = fESDpid->GetTOFResponse().GetStartTime(fmom); // T0best time
412     t0trackSigma = fESDpid->GetTOFResponse().GetStartTimeRes(fmom); // T0best time
413     Double_t sigma[5];
414     for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++){
415       sigma[ipart] = fESDpid->GetTOFResponse().GetExpectedSigma(fmom, inttime[ipart], AliPID::ParticleMass(ipart));}
416     sigmael=sigma[0];
417     sigmamu=sigma[1];
418     sigmapi=sigma[2];
419     sigmaka=sigma[3];
420     sigmapr=sigma[4];
421     Int_t TOFLAB[3];
422     track->GetTOFLabel(TOFLAB);
423     TOFlabel0=TOFLAB[0];
424     TOFlabel1=TOFLAB[1];
425     TOFlabel2=TOFLAB[2];
426     //TPCSignal = track->GetTPCsignal();
427     //TPCSigmaPI=fESDpid->NumberOfSigmasTPC(track,AliPID::kPion);
428     //TPCSigmaKA=fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon);
429     //TPCSigmaPR=fESDpid->NumberOfSigmasTPC(track,AliPID::kProton);
430     fTreeTrack->Fill(); 
431     
432     /*
433       kTOFout = TOF matching
434       kTIME = good integrated time
435       100000 > track->GetTOFsignal() > 12000 = TOF time reasanble range
436       tracklength > 365 = should be greater than the TOF radius (370 cm)
437     */
438     
439     
440     
441   
442     
443     
444   } // end of track loop
445   
446   frunOld=frun;
447   
448   
449   // Post output data  
450   PostData(1, fTreeTrack);
451   //PostData(1, fTreeEv);
452   PostData(2, hNumEv);
453
454 }      
455
456
457 //________________________________________________________________________
458 void TOFSpectrappAnalysis::Terminate(Option_t *) 
459 {
460   // Draw result to the screen
461   // Called once at the end of the query
462   Printf("*** CONSTRUCTOR CALLED ****");
463   fTreeTrack = dynamic_cast<TTree*> (GetOutputData(1));
464   if (!fTreeTrack) {
465     Printf("ERROR: fTreeTrack not available");
466     return;   
467   } 
468   
469   system("touch ok.job");
470 }
471
472
473 //________________________________________________________________________
474 Bool_t TOFSpectrappAnalysis::SelectOnImpPar(AliESDtrack* t) {
475   //
476   // cut on transverse impact parameter
477   //
478   Float_t d0z0[2],covd0z0[3];
479   t->GetImpactParameters(d0z0,covd0z0);
480   Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
481   Float_t d0max = 7.*sigma;
482   //
483   Float_t sigmaZ = 0.0146+0.0070/TMath::Power(t->Pt(),1.114758);
484   if (t->Pt() > 1) sigmaZ = 0.0216;
485   Float_t d0maxZ = 5.*sigmaZ;
486   //
487   if(TMath::Abs(d0z0[0]) < d0max && d0z0[1] < d0maxZ) return kTRUE;
488   return kFALSE;
489 }