[u/mrichter/AliRoot.git] / TPC / AliTPCcalibCosmic.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#include "Riostream.h"
#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TList.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TFile.h"

#include "AliTPCseed.h"
#include "AliESDVertex.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliESDInputHandler.h"

#include "AliTracker.h"
#include "AliMagFMaps.h"

#include "AliLog.h"

#include "AliTPCcalibCosmic.h"

#include "TTreeStream.h"
#include "AliTPCTracklet.h"

ClassImp(AliTPCcalibCosmic)


AliTPCcalibCosmic::AliTPCcalibCosmic() 
  :AliTPCcalibBase(),
  fHistNTracks(0),
  fClusters(0),
  fModules(0),
  fHistPt(0),
  fPtResolution(0),
  fDeDx(0)
{  
  AliInfo("Defualt Constructor");  
}


AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 
  :AliTPCcalibBase(),
  fHistNTracks(0),
  fClusters(0),
  fModules(0),
  fHistPt(0),
  fPtResolution(0),
  fDeDx(0)
{  
  SetName(name);
  SetTitle(title);
  AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
  AliTracker::SetFieldMap(field, kTRUE);  
  fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
  fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
  fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
  fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);  
  fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
  fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
  BinLogX(fDeDx);
  AliInfo("Non Defualt Constructor");  
}

AliTPCcalibCosmic::~AliTPCcalibCosmic(){
  //
  //
  //
}


void AliTPCcalibCosmic::Process(AliESDEvent *event) {
    
  if (!event) {
    Printf("ERROR: ESD not available");
    return;
  }
  
  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
  if (!ESDfriend) {
   Printf("ERROR: ESDfriend not available");
   return;
  }

  printf("Hallo world: Im here\n");
  
  Int_t n=event->GetNumberOfTracks(); 
  fHistNTracks->Fill(n);
  
  //track loop
  for (Int_t i=0;i<n;++i) { 
   AliESDtrack *track = event->GetTrack(i); 
   fClusters->Fill(track->GetTPCNcls());
   
   AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
   
   AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
   TObject *calibObject;
   AliTPCseed *seed = 0;
   for (Int_t l=0;calibObject=friendTrack->GetCalibObject(l);++l) {
    if (seed=dynamic_cast<AliTPCseed*>(calibObject)) break;
   }
 
   if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));

    
  }
  
  // dE/dx,pt and ACORDE study --> studies which need the pair selection  
  if (n > 2) return;
  
  for (Int_t i=0;i<n;++i) {
    AliESDtrack *track1 = event->GetTrack(i);
     
    Double_t d1[3];
    track1->GetDirection(d1);
    
    for (Int_t j=i+1;j<n;++j) {
     AliESDtrack *track2 = event->GetTrack(j);   
     Double_t d2[3];
     track2->GetDirection(d2);
       
     if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
     
      /*___________________________________ Pt resolution ________________________________________*/
      if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
       Double_t res = (track1->Pt() - track2->Pt());
       res = res/(2*(track1->Pt() + track2->Pt()));
       fPtResolution->Fill(100*res);
      }
      
      /*_______________________________ Propagation to ACORDE ___________________________________*/
      const Double_t AcordePlane = 850.; //distance of the central Acorde detectors to the beam line at y =0
      const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
     
      if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {        
       Double_t r[3];
       track1->GetXYZ(r);
       Double_t x,z;
       z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
       x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
       
       if (x > roof) {
        x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
        z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
       }
       if (x < -roof) {
        x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
        z = z -  TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
       }
       
       fModules->Fill(z, x);
      }
      
      if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
       Double_t r[3];
       track2->GetXYZ(r);
       Double_t x,z;
       z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
       x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
       
       if (x > roof) {
        x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
        z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));  
       }
       if (x < -roof) {
        x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
	z = z -  TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
       }       
       
       fModules->Fill(z, x);
      }
      
      printf("My stream level=%d\n",fStreamLevel);

      if (fStreamLevel>0){
	TTreeSRedirector * cstream =  GetDebugStreamer();
	printf("My stream=%p\n",cstream);
	if (cstream) {
	  (*cstream) << "Track" <<
	    "Track1.=" << track1 <<      // original track 1
	    "Track2.=" << track2 <<      // original track2
	    "\n";
	}
      }
  //     AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
//       AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
//       delete trackOut;
      

      break;            
     }     
    }
  }
  
  
}    


Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {

}


void AliTPCcalibCosmic::BinLogX(TH1 *h) {

  // Method for the correct logarithmic binning of histograms

  TAxis *axis = h->GetXaxis();
  int bins = axis->GetNbins();

  Double_t from = axis->GetXmin();
  Double_t to = axis->GetXmax();
  Double_t *new_bins = new Double_t[bins + 1];
   
  new_bins[0] = from;
  Double_t factor = pow(to/from, 1./bins);
  
  for (int i = 1; i <= bins; i++) {
   new_bins[i] = factor * new_bins[i-1];
  }
  axis->Set(bins, new_bins);
  delete new_bins;
  
}
Commit	Line	Data
f7f33dec	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	#include "Riostream.h"
	17	#include "TChain.h"
	18	#include "TTree.h"
	19	#include "TH1F.h"
	20	#include "TH2F.h"
	21	#include "TList.h"
	22	#include "TMath.h"
	23	#include "TCanvas.h"
	24	#include "TFile.h"
	25
	26	#include "AliTPCseed.h"
	27	#include "AliESDVertex.h"
	28	#include "AliESDEvent.h"
	29	#include "AliESDfriend.h"
	30	#include "AliESDInputHandler.h"
	31
	32	#include "AliTracker.h"
	33	#include "AliMagFMaps.h"
	34
	35	#include "AliLog.h"
	36
	37	#include "AliTPCcalibCosmic.h"
	38
	39	#include "TTreeStream.h"
	40	#include "AliTPCTracklet.h"
	41
	42	ClassImp(AliTPCcalibCosmic)
	43
	44
	45	AliTPCcalibCosmic::AliTPCcalibCosmic()
	46	:AliTPCcalibBase(),
	47	fHistNTracks(0),
	48	fClusters(0),
	49	fModules(0),
	50	fHistPt(0),
	51	fPtResolution(0),
	52	fDeDx(0)
	53	{
	54	AliInfo("Defualt Constructor");
	55	}
	56
	57
	58	AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t name, const Text_t title)
	59	:AliTPCcalibBase(),
	60	fHistNTracks(0),
	61	fClusters(0),
	62	fModules(0),
	63	fHistPt(0),
	64	fPtResolution(0),
65	fDeDx(0)
66	{
67	SetName(name);
68	SetTitle(title);
69	AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
70	AliTracker::SetFieldMap(field, kTRUE);
71	fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
72	fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
73	fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
74	fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);
75	fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
76	fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
77	BinLogX(fDeDx);
78	AliInfo("Non Defualt Constructor");
79	}
80
81	AliTPCcalibCosmic::~AliTPCcalibCosmic(){
82	//
83	//
84	//
85	}
86
87
88	void AliTPCcalibCosmic::Process(AliESDEvent *event) {
89
90	if (!event) {
91	Printf("ERROR: ESD not available");
92	return;
93	}
94
95	AliESDfriend ESDfriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
96	if (!ESDfriend) {
97	Printf("ERROR: ESDfriend not available");
98	return;
99	}
100
101	printf("Hallo world: Im here\n");
102
103	Int_t n=event->GetNumberOfTracks();
104	fHistNTracks->Fill(n);
105
106	//track loop
107	for (Int_t i=0;i<n;++i) {
108	AliESDtrack *track = event->GetTrack(i);
109	fClusters->Fill(track->GetTPCNcls());
110
111	AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
112
113	AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
114	TObject *calibObject;
115	AliTPCseed *seed = 0;
116	for (Int_t l=0;calibObject=friendTrack->GetCalibObject(l);++l) {
117	if (seed=dynamic_cast<AliTPCseed*>(calibObject)) break;
118	}
119
120	if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
121
122
123	}
124
125	// dE/dx,pt and ACORDE study --> studies which need the pair selection
126	if (n > 2) return;
127
128	for (Int_t i=0;i<n;++i) {
129	AliESDtrack *track1 = event->GetTrack(i);
130
131	Double_t d1[3];
132	track1->GetDirection(d1);
133
134	for (Int_t j=i+1;j<n;++j) {
135	AliESDtrack *track2 = event->GetTrack(j);
136	Double_t d2[3];
137	track2->GetDirection(d2);
138
139	if (d1[0]d2[0] + d1[1]d2[1] + d1[2]*d2[2] < -0.999) {
140
141	/___________________________________ Pt resolution ________________________________________/
142	if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
143	Double_t res = (track1->Pt() - track2->Pt());
144	res = res/(2*(track1->Pt() + track2->Pt()));
145	fPtResolution->Fill(100*res);
146	}
147
148	/_______________________________ Propagation to ACORDE ___________________________________/
149	const Double_t AcordePlane = 850.; //distance of the central Acorde detectors to the beam line at y =0
150	const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
151
152	if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {
153	Double_t r[3];
154	track1->GetXYZ(r);
155	Double_t x,z;
156	z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
157	x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
158
159	if (x > roof) {
160	x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
161	z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
162	}
163	if (x < -roof) {
164	x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
165	z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
166	}
167
168	fModules->Fill(z, x);
169	}
170
171	if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
172	Double_t r[3];
173	track2->GetXYZ(r);
174	Double_t x,z;
175	z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
176	x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
177
178	if (x > roof) {
179	x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
180	z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
181	}
182	if (x < -roof) {
183	x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
184	z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
185	}
186
187	fModules->Fill(z, x);
188	}
189
190	printf("My stream level=%d\n",fStreamLevel);
191
192	if (fStreamLevel>0){
193	TTreeSRedirector * cstream = GetDebugStreamer();
194	printf("My stream=%p\n",cstream);
195	if (cstream) {
196	(*cstream) << "Track" <<
197	"Track1.=" << track1 << // original track 1
198	"Track2.=" << track2 << // original track2
199	"\n";
200	}
201	}
202	// AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
203	// AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
204	// delete trackOut;
205
206
207
208
209
210	break;
211	}
212	}
213	}
214
215
216
217
218	}
219
220
221	Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
222
223	}
224
225
226	void AliTPCcalibCosmic::BinLogX(TH1 *h) {
227
228	// Method for the correct logarithmic binning of histograms
229
230	TAxis *axis = h->GetXaxis();
231	int bins = axis->GetNbins();
232
233	Double_t from = axis->GetXmin();
234	Double_t to = axis->GetXmax();
235	Double_t *new_bins = new Double_t[bins + 1];
236
237	new_bins[0] = from;
238	Double_t factor = pow(to/from, 1./bins);
239
240	for (int i = 1; i <= bins; i++) {
241	new_bins[i] = factor * new_bins[i-1];
242	}
243	axis->Set(bins, new_bins);
244	delete new_bins;
245
246	}
247