[u/mrichter/AliRoot.git] / TRD / qaRec / AliTRDmultiplicity.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.                  *
 **************************************************************************/


////////////////////////////////////////////////////////////////////////////
//                                                                        //
//  TRD multiplicity                                                      //
//                                                                        //
// Task to select true single tracks
// 
// Program flow:
// Part A - task AliTRDmultiplicity (within the framework qaRec):
//   For TRD standalone or TRD&TOF tracks I make a FitRiemanTilt
//   The resulting parameterization is dumped into a DebugStreamer written to a file
// 
// Part B – $ALICE_ROOT/TRD/qaRec/macros/TrackletsinTRD.C[h] (analysis macro):
//   The clusters are read in
//   The above produced file is read in
//   I make a fit through the parameterization of the  FitRiemanTilt -> in order to get a straight line (representing the track)
//   the distance of each cluster with respect to the straight line is calculated
//   If there is no cluster at a distance of 0.6  to 2 with respect to the track, the track is defined as a good track, i.e. clean single track
//                                                                        //
//  Authors: Yvonne C. Pachmayer <pachmay@physi.uni-heidelberg.de>        //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include <TClonesArray.h>
#include <TObjArray.h>
#include <TProfile.h>
#include <TH1.h>
#include <TH1F.h>
#include <TFile.h>
#include "AliTRDtracker.h"
#include "TTreeStream.h"

#include "AliPID.h"
#include "AliESDtrack.h"
#include "AliTrackReference.h"
#include "AliExternalTrackParam.h"
#include "AliTracker.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "AliTRDtrackerV1.h"
#include "AliTrackPointArray.h"
//#include "AliMagFMaps.h"
#include "AliTRDcluster.h"
#include "AliAnalysisManager.h"

#include "Cal/AliTRDCalPID.h"
#include "AliTRDgeometry.h"
#include "AliTRDmultiplicity.h"
#include "info/AliTRDtrackInfo.h"

ClassImp(AliTRDmultiplicity)

//____________________________________________________________________
AliTRDmultiplicity::AliTRDmultiplicity()
    :AliTRDrecoTask("Multiplicity", "Barrel Tracking Multiplicity")
  ,fEventCounter(0)
{
  //
  // Default constructor
    //
}

//____________________________________________________________________
AliTRDmultiplicity::~AliTRDmultiplicity()
{
}

//____________________________________________________________________
void  AliTRDmultiplicity::CreateOutputObjects()
{
  //
  // Create output objects
  //

  OpenFile(0, "RECREATE");

  TH1 *h = 0x0;
  fContainer = new TObjArray();
  for(Int_t is=0; is<AliTRDgeometry::kNsector; is++){
  fContainer->Add(h = new TH1F(Form("h_sector_%i", is), " ", 100,-10,10));
  }
} 

//____________________________________________________________________
void AliTRDmultiplicity::Exec(Option_t *)
{
  //
  // Do it
  //


  ULong_t status;
  AliTRDtrackInfo     *track = 0x0;
  AliExternalTrackParam *esd = 0x0;
  AliTRDtrackV1 *TRDtrack = 0x0;
  Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
  Int_t stand_alone=1;
  fEventCounter++;


  for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
    track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
    if(!track->HasESDtrack()) continue;
    status = track->GetStatus();

    AliTrackPoint points[6];
    Float_t xyz[3];
    memset(xyz, 0, sizeof(Float_t) * 3);
    Float_t tracklet_x[6];
    Float_t tracklet_y[6];
    Float_t tracklet_z[6];
    for(Int_t a=0;a<6;a++)
    {
        xyz[0] = x_anode[a];
        points[a].SetXYZ(xyz);
        tracklet_x[a]=-1000;
        tracklet_y[a]=-1000;
        tracklet_z[a]=-1000;
    }
    Int_t det_tracklet=600;


    // TRD standalone
    if(((status&AliESDtrack::kTRDout)>0) && !((status&AliESDtrack::kTRDin)>0))
    {
        TRDtrack = track->GetTrack();

        if(TRDtrack)
        {
            AliTRDseedV1 *tracklet = 0x0;
            Int_t counter_tracklet=0;
            for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++)
            {
                tracklet = TRDtrack->GetTracklet(itl);
                if(tracklet)
                {
                    if(tracklet->IsOK())
                    {
                        counter_tracklet++;
                        det_tracklet=tracklet->GetDetector();
                        //printf("%i %f %f \n",itl,tracklet->GetYfit(0),tracklet->GetZfit(0));
                        tracklet_x[itl]=tracklet->GetX0();
                        tracklet_y[itl]=tracklet->GetYfit(0);
                        tracklet_z[itl]=tracklet->GetZfit(0);
                    }
                }
            }
            // this cut is needed otherwise the order of tracklets in the fit function is not correct
            if(counter_tracklet==AliTRDgeometry::kNlayer) AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(TRDtrack), 0x0, kTRUE, counter_tracklet, points);
            else continue;

            if(fDebugLevel>=1){
              for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
                //printf("---------------------------------------- %i %i %f %f %f \n",counter_tracklet,il,points[il].GetX(),points[il].GetY(),points[il].GetZ());
                Double_t point_x=points[il].GetX();
                Double_t point_y=points[il].GetY();
                Double_t point_z=points[il].GetZ();


                (*fDebugStream)   << "TrackletsinTRD"
                    << "standalone=" << stand_alone
                    << "eventcounter=" << fEventCounter
                    << "layer="     << il
                    << "dettracklet=" << det_tracklet
                    << "xtracklet=" << tracklet_x[il]
                    << "xtrack="    << point_x
                    << "ytracklet=" << tracklet_y[il]
                    << "ytrack="    << point_y
                    << "ztracklet=" << tracklet_z[il]
                    << "ztrack="    << point_z
                    << "num_tracklets=" << counter_tracklet
                    << "\n";

              }
            }
        }
    } // end TRD standalone


    // TPC cluster selection for cosmic data
    if((track->GetTPCncls())<40) continue;
    // missing TPC propagation
    if(!(status&AliESDtrack::kTPCout)) continue;
    stand_alone=0;

    esd = track->GetESDinfo()->GetOuterParam();
    // calculate sector - does not matter if track ends in TPC or TRD - sector number independent
    Double_t alpha;
    alpha=-100;
    Int_t fSector=100;
    if(esd){
        alpha=esd->GetAlpha();
        fSector = (Int_t)((alpha+ TMath::Pi())/(20*TMath::Pi()/180));
        if((fSector>-1) && (fSector<18))
        {
        }
        else
        {
            continue;
        }
    }


    // only these sectors have a TRD detector at the moment
   if(fSector==0||fSector==8||fSector==9||fSector==17)
    {
        TRDtrack = track->GetTrack();
        if(!TRDtrack) continue;
        AliTRDseedV1 *tracklet = 0x0;
        Int_t counter_tracklet=0;


        for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
            tracklet = TRDtrack->GetTracklet(itl);
            if(!tracklet || !tracklet->IsOK()) continue;
            counter_tracklet++;
            det_tracklet=tracklet->GetDetector();
            tracklet_x[itl]=tracklet->GetX0();
            tracklet_y[itl]=tracklet->GetYfit(0);
            tracklet_z[itl]=tracklet->GetZfit(0);

            /*
            AliTRDcluster *c = 0x0;
            for(Int_t itime = 0; itime < AliTRDseedV1::kNtb; itime++){
                c = tracklet->GetClusters(itime);
                if(!c) continue;
//                Float_t cluster_x = c->GetX();
//                Bool_t isprop;
//                isprop=kFALSE;
//                isprop=AliTracker::PropagateTrackTo(esd,(Double_t)cluster_x,0.105,5,kFALSE);
//                if(isprop)
//                {
//                    Int_t detector = c->GetDetector();
//                    ((TH1F*)fContainer->At(fSector))->Fill((c->GetY())-(esd->GetY()));
//                    if(c->GetY()-esd->GetY()>);
//                    printf("diff: %f\n",c->GetY()-(esd->GetY()));
//                }
            }
            */

        } // loop over tracklets ends

        if(counter_tracklet==AliTRDgeometry::kNlayer)
        {
            AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(TRDtrack), 0x0, kTRUE, counter_tracklet, points);
        }
        else continue;

       
        if(fDebugLevel>=1){
          for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
            //printf("---------------------------------------- %i %i %f %f %f \n",counter_tracklet,il,points[il].GetX(),points[il].GetY(),points[il].GetZ());
            Double_t point_x=points[il].GetX();
            Double_t point_y=points[il].GetY();
            Double_t point_z=points[il].GetZ();
            
            (*fDebugStream)   << "TrackletsinTRD"
                << "standalone=" << stand_alone
                << "eventcounter=" << fEventCounter
                << "layer="     << il
                << "dettracklet=" << det_tracklet
                << "xtracklet=" << tracklet_x[il]
                << "xtrack="    << point_x
                << "ytracklet=" << tracklet_y[il]
                << "ytrack="    << point_y
                << "ztracklet=" << tracklet_z[il]
                << "ztrack="    << point_z
                << "num_tracklets=" << counter_tracklet
                << "\n";
          }
        }

    }
  }

  
  PostData(0, fContainer);
}

//____________________________________________________________________
void AliTRDmultiplicity::Terminate(Option_t *)
{
  //
  // Terminate
  //

  if(fDebugStream){ 
    delete fDebugStream;
    fDebugStream = 0x0;
    fDebugLevel = 0;
  }

  fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
  if (!fContainer) {
    Printf("ERROR: list not available");
    return;
  }

}
Commit	Line	Data
c792ca2e	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	// TRD multiplicity //
	20	// //
	21	// Task to select true single tracks
	22	//
	23	// Program flow:
	24	// Part A - task AliTRDmultiplicity (within the framework qaRec):
	25	// For TRD standalone or TRD&TOF tracks I make a FitRiemanTilt
	26	// The resulting parameterization is dumped into a DebugStreamer written to a file
	27	//
	28	// Part B – $ALICE_ROOT/TRD/qaRec/macros/TrackletsinTRD.C[h] (analysis macro):
	29	// The clusters are read in
	30	// The above produced file is read in
	31	// I make a fit through the parameterization of the FitRiemanTilt -> in order to get a straight line (representing the track)
	32	// the distance of each cluster with respect to the straight line is calculated
	33	// If there is no cluster at a distance of 0.6 to 2 with respect to the track, the track is defined as a good track, i.e. clean single track
	34	// //
	35	// Authors: Yvonne C. Pachmayer <pachmay@physi.uni-heidelberg.de> //
	36	// //
	37	////////////////////////////////////////////////////////////////////////////
	38
	39	#include <TClonesArray.h>
	40	#include <TObjArray.h>
	41	#include <TProfile.h>
	42	#include <TH1.h>
	43	#include <TH1F.h>
	44	#include <TFile.h>
	45	#include "AliTRDtracker.h"
	46	#include "TTreeStream.h"
	47
	48	#include "AliPID.h"
	49	#include "AliESDtrack.h"
	50	#include "AliTrackReference.h"
	51	#include "AliExternalTrackParam.h"
	52	#include "AliTracker.h"
	53	#include "AliTRDseedV1.h"
	54	#include "AliTRDtrackV1.h"
	55	#include "AliTRDtrackerV1.h"
	56	#include "AliTrackPointArray.h"
	57	//#include "AliMagFMaps.h"
	58	#include "AliTRDcluster.h"
	59	#include "AliAnalysisManager.h"
	60
	61	#include "Cal/AliTRDCalPID.h"
	62	#include "AliTRDgeometry.h"
	63	#include "AliTRDmultiplicity.h"
873458ab	64	#include "info/AliTRDtrackInfo.h"
c792ca2e	65
	66	ClassImp(AliTRDmultiplicity)
	67
	68	//____________________________________________________________________
	69	AliTRDmultiplicity::AliTRDmultiplicity()
	70	:AliTRDrecoTask("Multiplicity", "Barrel Tracking Multiplicity")
	71	,fEventCounter(0)
	72	{
	73	//
	74	// Default constructor
	75	//
	76	}
	77
	78	//____________________________________________________________________
	79	AliTRDmultiplicity::~AliTRDmultiplicity()
	80	{
	81	}
	82
	83	//____________________________________________________________________
	84	void AliTRDmultiplicity::CreateOutputObjects()
	85	{
	86	//
	87	// Create output objects
	88	//
	89
	90	OpenFile(0, "RECREATE");
	91
	92	TH1 *h = 0x0;
	93	fContainer = new TObjArray();
	94	for(Int_t is=0; is<AliTRDgeometry::kNsector; is++){
	95	fContainer->Add(h = new TH1F(Form("h_sector_%i", is), " ", 100,-10,10));
	96	}
	97	}
	98
	99	//____________________________________________________________________
	100	void AliTRDmultiplicity::Exec(Option_t *)
	101	{
	102	//
	103	// Do it
	104	//
	105
	106
	107
	108	ULong_t status;
	109	AliTRDtrackInfo *track = 0x0;
	110	AliExternalTrackParam *esd = 0x0;
	111	AliTRDtrackV1 *TRDtrack = 0x0;
	112	Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
	113	Int_t stand_alone=1;
	114	fEventCounter++;
	115
	116
	117	for(Int_t itrk=0; itrk<fTracks->GetEntriesFast(); itrk++){
	118	track = (AliTRDtrackInfo*)fTracks->UncheckedAt(itrk);
	119	if(!track->HasESDtrack()) continue;
	120	status = track->GetStatus();
	121
	122	AliTrackPoint points[6];
	123	Float_t xyz[3];
	124	memset(xyz, 0, sizeof(Float_t) * 3);
	125	Float_t tracklet_x[6];
	126	Float_t tracklet_y[6];
	127	Float_t tracklet_z[6];
	128	for(Int_t a=0;a<6;a++)
129	{
130	xyz[0] = x_anode[a];
131	points[a].SetXYZ(xyz);
132	tracklet_x[a]=-1000;
133	tracklet_y[a]=-1000;
134	tracklet_z[a]=-1000;
135	}
136	Int_t det_tracklet=600;
137
138
139
140	// TRD standalone
141	if(((status&AliESDtrack::kTRDout)>0) && !((status&AliESDtrack::kTRDin)>0))
142	{
143	TRDtrack = track->GetTrack();
144
145	if(TRDtrack)
146	{
147	AliTRDseedV1 *tracklet = 0x0;
148	Int_t counter_tracklet=0;
149	for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++)
150	{
151	tracklet = TRDtrack->GetTracklet(itl);
152	if(tracklet)
153	{
154	if(tracklet->IsOK())
155	{
156	counter_tracklet++;
157	det_tracklet=tracklet->GetDetector();
158	//printf("%i %f %f \n",itl,tracklet->GetYfit(0),tracklet->GetZfit(0));
159	tracklet_x[itl]=tracklet->GetX0();
160	tracklet_y[itl]=tracklet->GetYfit(0);
161	tracklet_z[itl]=tracklet->GetZfit(0);
162	}
163	}
164	}
165	// this cut is needed otherwise the order of tracklets in the fit function is not correct
166	if(counter_tracklet==AliTRDgeometry::kNlayer) AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(TRDtrack), 0x0, kTRUE, counter_tracklet, points);
167	else continue;
168
169	if(fDebugLevel>=1){
170	for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
171	//printf("---------------------------------------- %i %i %f %f %f \n",counter_tracklet,il,points[il].GetX(),points[il].GetY(),points[il].GetZ());
172	Double_t point_x=points[il].GetX();
173	Double_t point_y=points[il].GetY();
174	Double_t point_z=points[il].GetZ();
175
176
177	(*fDebugStream) << "TrackletsinTRD"
178	<< "standalone=" << stand_alone
179	<< "eventcounter=" << fEventCounter
180	<< "layer=" << il
181	<< "dettracklet=" << det_tracklet
182	<< "xtracklet=" << tracklet_x[il]
183	<< "xtrack=" << point_x
184	<< "ytracklet=" << tracklet_y[il]
185	<< "ytrack=" << point_y
186	<< "ztracklet=" << tracklet_z[il]
187	<< "ztrack=" << point_z
188	<< "num_tracklets=" << counter_tracklet
189	<< "\n";
190
191	}
192	}
193	}
194	} // end TRD standalone
195
196
197
198	// TPC cluster selection for cosmic data
199	if((track->GetTPCncls())<40) continue;
200	// missing TPC propagation
201	if(!(status&AliESDtrack::kTPCout)) continue;
202	stand_alone=0;
203
204	esd = track->GetESDinfo()->GetOuterParam();
205	// calculate sector - does not matter if track ends in TPC or TRD - sector number independent
206	Double_t alpha;
207	alpha=-100;
208	Int_t fSector=100;
209	if(esd){
210	alpha=esd->GetAlpha();
211	fSector = (Int_t)((alpha+ TMath::Pi())/(20*TMath::Pi()/180));
212	if((fSector>-1) && (fSector<18))
213	{
214	}
215	else
216	{
217	continue;
218	}
219	}
220
221
222
223	// only these sectors have a TRD detector at the moment
224	if(fSector==0\|\|fSector==8\|\|fSector==9\|\|fSector==17)
225	{
226	TRDtrack = track->GetTrack();
227	if(!TRDtrack) continue;
228	AliTRDseedV1 *tracklet = 0x0;
229	Int_t counter_tracklet=0;
230
231
232	for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
233	tracklet = TRDtrack->GetTracklet(itl);
234	if(!tracklet \|\| !tracklet->IsOK()) continue;
235	counter_tracklet++;
236	det_tracklet=tracklet->GetDetector();
237	tracklet_x[itl]=tracklet->GetX0();
238	tracklet_y[itl]=tracklet->GetYfit(0);
239	tracklet_z[itl]=tracklet->GetZfit(0);
240
241	/*
242	AliTRDcluster *c = 0x0;
243	for(Int_t itime = 0; itime < AliTRDseedV1::kNtb; itime++){
244	c = tracklet->GetClusters(itime);
245	if(!c) continue;
246	// Float_t cluster_x = c->GetX();
247	// Bool_t isprop;
248	// isprop=kFALSE;
249	// isprop=AliTracker::PropagateTrackTo(esd,(Double_t)cluster_x,0.105,5,kFALSE);
250	// if(isprop)
251	// {
252	// Int_t detector = c->GetDetector();
253	// ((TH1F*)fContainer->At(fSector))->Fill((c->GetY())-(esd->GetY()));
254	// if(c->GetY()-esd->GetY()>);
255	// printf("diff: %f\n",c->GetY()-(esd->GetY()));
256	// }
257	}
258	*/
259
260	} // loop over tracklets ends
261
262	if(counter_tracklet==AliTRDgeometry::kNlayer)
263	{
264	AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(TRDtrack), 0x0, kTRUE, counter_tracklet, points);
265	}
266	else continue;
267
268
269
270	if(fDebugLevel>=1){
271	for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
272	//printf("---------------------------------------- %i %i %f %f %f \n",counter_tracklet,il,points[il].GetX(),points[il].GetY(),points[il].GetZ());
273	Double_t point_x=points[il].GetX();
274	Double_t point_y=points[il].GetY();
275	Double_t point_z=points[il].GetZ();
276
277	(*fDebugStream) << "TrackletsinTRD"
278	<< "standalone=" << stand_alone
279	<< "eventcounter=" << fEventCounter
280	<< "layer=" << il
281	<< "dettracklet=" << det_tracklet
282	<< "xtracklet=" << tracklet_x[il]
283	<< "xtrack=" << point_x
284	<< "ytracklet=" << tracklet_y[il]
285	<< "ytrack=" << point_y
286	<< "ztracklet=" << tracklet_z[il]
287	<< "ztrack=" << point_z
288	<< "num_tracklets=" << counter_tracklet
289	<< "\n";
290	}
291	}
292
293	}
294	}
295
296
297	PostData(0, fContainer);
298	}
299
300	//____________________________________________________________________
301	void AliTRDmultiplicity::Terminate(Option_t *)
302	{
303	//
304	// Terminate
305	//
306
307	if(fDebugStream){
308	delete fDebugStream;
309	fDebugStream = 0x0;
310	fDebugLevel = 0;
311	}
312
313	fContainer = dynamic_cast<TObjArray*>(GetOutputData(0));
314	if (!fContainer) {
315	Printf("ERROR: list not available");
316	return;
317	}
318
319	}
320
321
322